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The problem of the determination of the nuclear surface and surface symmetry energy is addressed in the 
framework of the Extended Thomas Fermi (ETF) approximation using Skyrme functionals. We propose an 
analytical model for the density profiles with variationally determined diffuseness parameters. For the case 
of symmetric nuclei, the resulting ETF functional can be exactly integrated, leading to an analytical formula 
expressing the surface energy as a function of the couplings of the energy functional. The importance of non¬ 
local terms is stressed, which cannot be simply deduced from the local part of the functional. In the case of 
asymmetric nuclei, we propose an approximate expression for the diffuseness and the surface energy. These 
quantities are analytically related to the parameters of the energy functional. In particular, the influence of the 
different equation of state parameters can be explicitly quantified. Detailed analyses of the different energy 
components (local/non-local, isoscalar/isovector, surface/curvature and higher order) are also performed. Our 
analytical solution of the ETF integral improves over previous models and leads to a precision better than 
200 keV per nucleon in the determination of the nuclear binding energy for dripline nuclei. 


I. INTRODUCTION 

Skyrme functionals have been widely used to describe nu¬ 
clear structure properties, with different level of sophistica¬ 
tion in the many-body treatment, from the simplest Thomas- 
Fermi cii to modern multi-reference calculations |2|. The 
most basic observable accessible to the functional treatment 
is given by nuclear mass, allowing the analysis of the differ¬ 
ent mass components in terms of bulk and surface properties, 
as well as isovector and isoscalar properties. The theoretical 
prediction of nuclear mass is not only important in itself, but it 
is also a fundamental tool to optimize the different functional 
forms and associated parameters, for an increasing predictive 
power of density functional calculations H]. Indeed mass pre¬ 
dictions from microcopic density functionals nowadays starts 
to equalize the most precise phenomenological mass formulas 
available in the 1 iteratiire14 6|. 

For practical applications in nuclear structure or nuclear 
astrophysics problems, different parametrizations of nuclear 
masses fitted on density functional calculations with Skyrme 
forces have been proposed The limitation of these 

works is that the different coefficients are not analytically cal¬ 
culated but they result from the fit to the numerically deter¬ 
mined nuclear masses. As a consequence, the fit has to be 
performed again each time that the functional is improved by 
adding further constraints from the rapidly improving experi¬ 
mental data. Moreover, the absence of an analytical link be¬ 
tween the Skyrme parameters and the coefficients of the mass 
formula implies that it is difficult to make an unambiguous 
correlation between the different parts of the mass functional 
and the physical properties of the effective interaction. For 
these reasons, it appears interesting to search for an analytical 
expression of the mass formula coefficients, directly linked to 
the functional form and parameters of the Skyrme interaction. 
The derivation of such an analytical formula is the purpose of 
this paper. 

An especially appealing formalism when seeking for an¬ 
alytical expressions is the semi-classical Extended-Thomas- 


Fermi (ETF) approach, which is based on an expansion in 
powers of h of the energy functional 19.!, HH-EH] • The advan¬ 
tage of the ETF approximation is that the non-local terms in 
the energy density functional are entirely replaced by local 
gradients. As a consequence, the energy functional solely de¬ 
pends on the local particle densities. Thus, the energy of any 
arbitrary nuclear configuration can be calculated if the neu¬ 
tron and proton density profiles p„ and p p are given through 
a parametrized form. These density profiles are those of the 
ground state, or of any arbitrary excited state. A large number 
of configurations can therefore be explored, and this appealing 
property of ETF has been used to study nuclear configurations 
in dilute stellar matter contributing to the sub-saturation finite 
temperature equation of state II17I - U9I1 . On the other side, the 
well known limitation of ETF is that only the smooth part of 
the nuclear mass can be addressed, and shell effects have to be 
added on top, for instance through the well known Strutinsky 
integral theoremfTlj] ■ 

In this paper, we will consider an ETF expansion up to the 
second h 2 -order, and limit ourselves to the smooth part of the 
mass functional. 

The plan of the paper is as follows. 

Section m addresses the problem of symmetric nuclei. A 
single density profile is supposed for protons and neutrons, 
and symmetry breaking effects are included by accounting 
for the Coulomb modification of the bulk compressibility. In 
this simplified case, the ETF integrals can be analytically in¬ 
tegrated leading to a very transparent form for the surface and 
curvature terms of the nuclear energy (section III Al l and of 
the surface diffuseness (section Hi Bl) . In this same section we 
also retrive (section Hi Cl) that in a one-dimensional geometry 
the local and non-local terms are related, and the surface ten¬ 
sion can be consequently be expressed as a function of the 
local terms only [12]. This means that the surface tension 
solely depends on the local components of the energy density 
functional, that is the bulk properties of nuclear matter, and 
it does not depend on the non-local gradient and spin-orbit 
terms. This remarkable property however breaks down in 
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spherical symmetry, and any, even slight approximation to the 
exact variational profile, for instance the use of parametrized 
densities, increases the difference between local and non-local 
contributions to the surface energy. As a consequence, using 
parametrized density profiles, the contribution of non-local 
terms has to be carefully calculated independently of the local 
part, and the two separate contributions must be summed up 
to obtain the surface energy and the surface tension. 

The more general problem of isospin asymmetric nuclei is 
studied in section[III] We first demonstrate in section llll Al that 
a large part of the isospin dependence can be accounted for, 
if the asymmetry dependence of the saturation density is in¬ 
troduced in the nuclear bulk. The residual surface symmetry 
part is then defined in terms of the isovector density. This en¬ 
ergy density term is not analytically integrable, meaning that 
approximations have to be performed. We propose in sec¬ 
tion UnB] two different approximations and critically discuss 
their validity in comparison both to numerical integration of 
the ETF functional, and to complete Hartree-Fock (HF) cal¬ 
culations using the same functional (section llll B 3b . The first 
approximation, inspired from Ref. lllql . consists in neglect¬ 
ing the neutron skin (section flllB II) . Surprisingly enough, 
this very crude approximation leads to an overestimation of 
Hartree-Fock energies of medium-heavy and heavy nuclei of 
no more than 200-400 keV/nucleon even for the most extreme 
dripline nuclei. Again, such an accuracy can be obtained only 
if both local and non-local terms in the energy functional are 
separately calculated, meaning that the symmetry energy does 
not only depend on bulk nuclear matter properties. This might 
be at the origin of the well known ambiguities in the extraction 
of the symmetry energy from density functional calculations 
of finite nuclear properties flOl I36L142114311 . A better accuracy 
for neutron rich nuclei is obtained if isospin fluctuations are 
accounted for, and in section llll B 2l the approximation is made 
that the surface symmetry energy density is strongly peaked at 
the nuclear surface. 

Finally the complete mass formula is calculated for differ¬ 
ent representative Skyrme functionals in section[IV] The qual¬ 
itative behavior of the different energy components, that is the 
surface, curvature and higher order terms decomposed into 
isovector and isoscalar parts, and local and non-local parts, 
is discussed. The different analytical expressions for the mass 
functional are explicitly demonstrated in the appendix and can 
be readily used with any Skyrme interaction. The paper is 
summarized in section[V] and conclusions are given. 


II. SYMMETRIC NUCLEI 


Let us first consider a locally symmetric matter distribu¬ 
tion, that is characterized by a single density profile which is 
supposed to be identical for protons and neutrons. This ide¬ 
alized situation is not completely realistic even in N = Z nu¬ 
clei because of isospin symmetry breaking due to Coulomb. 
However, it was shown |[19;] that a great part of the Coulomb 
isospin symmetry breaking effects can be included simply ac¬ 
counting for the Coulomb modification of the bulk compress¬ 
ibility IfM [27ll35ll . This single-density model leads to an ex¬ 


cellent reproduction of the microscopically calculated as well 
as experimental!y_measured magic N = Z nuclei over the pe¬ 
riodic table fT9[ [ 2 OI. 

The idealized case of a common density profile for protons 
and neutrons has the advantage of leading to exact formu¬ 
las for the nuclear binding energy, as we explicitly show in 
this section. As we will see, this allows disentangling in a 
non-ambiguous way bulk, surface, curvature as well as higher 
order terms, and to determine exact relations connecting the 
different energy components to the parameters of the energy 
functional. 

Neglecting spin-gradient terms, the ETF Skyrme energy 
density reads, 


Jt?\p]=h{p)- 


2m* (p) 


t 2 - 


C /in -|^pm*(p))(Vp) 2 . 

(1) 


In this expression, m*{p) is the density dependent effective 
mass, m/m* = 1 + j/rC e ffp, the kinetic energy density con¬ 
sists of a zero order Thomas-Fermi term To as well as of a 
second order local and non-local correction t 2 = X l 2 + t^ 1 : 
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with / = m/m*. The local terms are given by: 

' ,{p) = i^-p l+c " p2+ClP °* 2 ’ (5) 


and gradient terms arise both from the non-local and the spin- 
orbit part of the Skyrme functional. (Co,C^,C e ff,Cfi n ,C so ,a) 
are Skyrme parameters, given in appendix [A] Spin-gradient 
terms are not considered in the applications of this paper, but 
their inclusion is straightforward. Full expressions are given 
in appendix [A] We will also limit ourselves to spherical sym¬ 
metry throughout the paper. 

To compute Eq. £[]), the density profile p(r) is required. 
The most common choice in the literature |2Q] consists in tak¬ 
ing a Fermi function. In particular, it was shown lfl9ll that a 
Fermi function succeeds in well reproducing the density pro¬ 
files and the corresponding energy calculated with the spheri¬ 
cal HF model. The density profile reads. 


P0) = p sa ,F (r) ; F(r) = (l W r “ S)/ ") 1 . (6) 

In this equation, p sat is the saturation density of symmetric 
nuclear matter, and R is the radius parameter related to the 
particle number of the nucleus 






(7) 


Let us observe that Eq. ([7} is a finite expansion and does not 
require any assumption except that e R i a <C 1, that is a < R 
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|‘22f|. If, in addition, we assume a we can invert equation 
O to get at the fourth order in (a/R) 



where Rhs = A 1 V vfl/ is the equivalent homogeneous sphere 

radius, and r sat = (ixpsat) is the mean radius per nu¬ 
cleon. 

The two other parameters entering Eq. © are the diffuse¬ 
ness a of the density profile, which is analytically derived in 
section lHBl and the saturation density p sat which corresponds 
to the equilibrium density of homogeneous infinite symmetric 
matter. 


A. Ground state energy 

Integrating in space Eq. 0} computed with the parametrized 
density profile given by Eq. ©, allows obtaining the total en¬ 
ergy £ of a nucleus of a mass A defined by Eq. £7]): 

E= f°° dr J4?[p(r)\. (9) 

Jo 

Within the nucleus total binding energy, it is interesting to dis¬ 
tinguish the bulk, surface, and curvature components corre¬ 
sponding to different functional dependences on the nuclear 
size l2lh . 

The bulk energy E b is the energy of a finite volume of nu¬ 
clear matter. It corresponds to the energy that the nucleus 
would have without finite-size effects, defined by: 

E b =J%atV HS = KatA , (10) 

where Vjjs = 4/3717?^ = A/p sat is the equivalent homoge¬ 
neous sphere volume and A s „, is the energy per particle at sat¬ 
uration: 

= Ckin — —Psm + Co Psat + C 3 Psat * i (11) 

Psat 

with Cki„ = |/r/(2m)(37r 2 /2) 2 / 3 and m* at = m*(p sat ). The 
finite-size correction to the bulk energy E s is defined as the 
total energy after the bulk is removed, that is 

E s = J dr Jt?\p{r)}-J% a ,V HS 

poo 

= 47 ij d r{jt°\p(r)]-\ sa ,p(r)}r 2 . (12) 

This finite size contribution E s will be called the surface en¬ 
ergy in the following. Let us notice that if we have properly 
removed the bulk energy part by Eq. ( flOb . the surface energy 
should scale with A with a dependence slower than linear, but 
the dependence can be different from A 2 / 3 because of curva¬ 
ture and higher order terms, see also ref .1I2TI1. We will see in 



the following that it is indeed the case in spherical symmetry. 
In the energy density Jff’lp] Eq. 0}. we can distinguish the 
non-local terms which depend on the density derivatives and 
are pure finite-size effects, from the local energy part hip) 
which only depends on the equation of state and on the density 
profile. We then write the surface energy as E s = E/ + E' v/ % 
with E^ the local part and E^ L the non-local one 1164(1 : 

E s= 4n f d rih[p{r)\- h ^ Psc "^ p(r)\r 2 , (13) 

Jo L Psat J 

(14) 


To obtain Eq. (fl4l >. we have changed the Laplace deriva¬ 
tives into gradients in the kinetic part, see Eqs. (1A16I) . in¬ 
tegrating by parts. Making a simple variable change, the 
originally 3-dimensional integral can be turned into the sum 
of three 1-dimensional integrals (see appendix [B]). Then a 
very accurate approximation, that is with an error less than 
(exp(— 5a/3R) — exp(— a/R)), allows to analytically inte¬ 
grate the differences of Fermi functions, such that the local 
and non-local terms can be written as a function of the effec¬ 
tive interaction parameters as (calculation details are given in 
appendixICl): 
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where the coefficients ,, depend on the sat- 

surf(curv)(ind ) “ 

uration density p sat and on the Skyrme parameters 
Co,C 3 ,C e ff,a,Cfj n ,C so , and where we have anticipated the 
(slight) A-dependence of the diffuseness in the most general 
case (see section HTbI i. 

The coefficients %’[' and corresponding to the local and 
non-local energy components read (see appendix IC It : 
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where 8m sat = (m — m* at )/m* at , V so = —mC 2 0 /ti 2 , and where 

(k) 

we have introduced the coefficients pj, defined by equa¬ 
tion (IB2b . Their numerical values are given in the same 
appendix. In order to have an analytical expression, we 
have made in eqs. (I20l )-(l22l) a Taylor expansion of the effec¬ 
tive mass inverse / _1 = (—l)'(5m)‘. This expansion is 

rapidly convergent: a truncation at i max = 7 produces an error 
~ 1% at the highest possible density p sat = 0.16 fm - in the 
case of the SLy4 interaction. 

Equations ( fiTt i and (fl6t show that the dominant surface ef¬ 
fect in the symmetric nucleus energetics is, as expected, a term 
oc A 2 / 3 . As it is well known, this term fully exhausts the finite- 
size effects given by the presence of a nuclear surface in the 
one-dimensional case of a semi-infinite slab geometry mu. 
Indeed in this case we have, see amtendix lB 21 

Ef ab = / + “d*{^[p(*)] - A saI p(x)}. (23) 

The evaluation of the integral Eq. (1231) leads to the same <=c 
A 2 / 3 term as in the spherical geometry, with a modified form 
factor 4-nRjjy 

a = a L + a NL = (tf s L urf + ■ (24) 

where a = lim A^ >00 Ef ab /A 2 ! 2 is the slab surface tension. The 
form factor difference between the surface energy of the slab 


and the one in spherical symmetry signs the difference of ge¬ 
ometry, and the spherical surface energy is the surface area 
multiplied by the energy per unit area of the infinite tangent 
plane. Let us notice that since the mass cannot be defined in 
the semi-infinite medium, the diffuseness in Eq. (l24l > is a con¬ 
stant. 

In a three-dimensional geometry, the existence of a surface 
leads to additional finite-size terms, even in the spherically 
symmetric case, as shown by eqs. (IT5l) . ifiTil i. The terms pro¬ 
portional to A 1 / 3 are the so-called curvature terms which cor¬ 
rect the surface energy with respect to the slab tangent limit. It 
is interesting to notice that we also have A-independent terms, 
which are rarely accounted for in the literature but turn out to 
be important for light nuclei m. As it can be seen in Eq. ([8]), 
higher order terms are of the order A -1 / 3 and are systemat¬ 
ically neglected in this work. This Taylor expansion is known 
in the literature as the leptodermous expansion l2ll [381 . It 
is interesting to observe that both local and non-local plane 
surface, curvature, and mass independent energy components 
arise even if no explicit gradient term is included in the func¬ 
tional. As a consequence, surface properties are determined 
by a complex interplay between equation of state properties 
and specific finite nuclei properties like spin-orbit and finite 
range. Using the definitions of the energy per particle at satu¬ 
ration X sa t = h/p | = dh/dp \ psai and the nuclear symmet¬ 
ric matter incompressibility K sat =9p 2 at d 2 (h/p)/dp 2 \ Psat , we 
can express the local energy eqs. <fl7b . (fill) and ( ITbl i as a func¬ 
tion of nuclear matter properties only, using the following ex¬ 
pressions: 
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The expression of the coefficients ^ N>L greatly simplifies 
if we consider a simplistic Zamick-type interaction, with a = 
1 and m = m*: 


^Lrf — (^5 ^ 5/3 + T^) e sat ~ Tj^sat, (28) 

^ u ‘f = 24 iTn + 2 CfinPsat ’ ’ (29) 

where we have introduced the Fermi energy per nucleon at 
saturation e F sat = § C kin p 2 /, 3 . 

We can see that even in this oversimplified model the nu¬ 
clear surface properties cannot be simply reduced to EoS pa¬ 
rameters. We can also gather the local and non-local terms 
in order to classify finite-size effects according to the rank of 
the Taylor expansion. Thus we introduce the surface E sur f, 
curvature E curv and ,4-indepcndcnt E m ,i energy components: 
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(30) 
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(32) 


We can see that all terms are multiplied by a power of the 
diffuseness except the non-local curvature part which is not. 
The role of the diffuseness on the surface properties thus de¬ 
pends on the rank of the Taylor expansion (surface, curva¬ 
ture, independent,...) and is not the same for the local or the 
non-local part. The functional difference between the local 
and non-local terms comes from the squared density gradi¬ 
ent appearing in the non-local energy Eq. (IT4l> . Globally, if 
the diffuseness is high, the local energy dominates over the 
non-local one, see eqs. (f30l i-([32li. This is easy to understand: 
in the limit of a purely local energy functional, the optimal 
configuration corresponds to a homogeneous hard sphere at 
saturation density, given by a = 0. The existence of a finite 
diffuseness for atomic nuclei is due to the presence of non¬ 
local terms in the functional, because of both explicit gradient 
interactions and of quantum effects on the kinetic energy den¬ 
sity. Let us notice that both effects are present even in the 
simplified eqs. d28l . d29l >. 


B. Analytical expression for the diffuseness 


The ground state energy of this model for symmetric nu¬ 
clei is given by the minimization of the energy per nucleon 
8(E/A) = 0 with the constraint of a given mass number A. 
We have seen in section iHAl that the only unconstrained pa¬ 
rameter of the model is the diffuseness parameter a. Though 
it does not play any role in the bulk energy, it is an essential 
ingredient for the surface energy E s given by eqs. o. eg. 
The diffuseness parameter can therefore be obtained from the 
variational equation |fl2| |: 


dE s 

da 


= 0 . 


(33) 


In principle, one should also add the surface Coulomb energy 
into E s , which would change the variational equation. How¬ 
ever, the resulting correction on a is very small & 

Equation (l33l) turns out to be particularly simple in the one¬ 
dimensional case of semi-infinite matter, or equivalently ne¬ 
glecting curvature and A-independent terms when consider¬ 
ing nuclei. Indeed, in this case, Eq. (l24l) leads to an analytical 
solution, already obtained in Ref. 11211 : 


a = 


cupNL 
®surf 



(34) 


This equation shows that the slab diffuseness a is determined 
by the balance between the local terms, which favour low dif¬ 
fuseness values corresponding to a hard sphere of matter at 
saturation density; and non-local terms which favour a large 
diffuseness corresponding to matter close to uniformity. 

The complete spherical case leads to the following 4 th de¬ 
gree polynomial equation: 

3 <d(—) +2 K L un A 1/3 (—) (35) 

\ T sat J \ ^sat J 

+ (*i+/3+4-o) (f) 2 - 4-0 a>/3=0. 

V “sat / \ 'sat J ^sat 

which has to be solved numerically. 

Notice that the coefficient c tff l F v does not contribute to this 
equation since, as already mentioned, it does not depend on a , 
cf. Eq. iOTI) . The solution of this equation, as well as the slab 
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FIG. 1. (Color online) Diffuseness (upper panel) and energy 
per nucleon (lower panel) of symmetric nuclei as a function of the 
mass number. Full red lines: calculations using the slab diffuseness 
Eq. {B Dashed blue lines: calculations using the spherical dif¬ 
fuseness Eq. 03 Dash-dotted green lines: calculations using the 
diffuseness fitted from HF density profiles (20). Star symbols: full 
Hartree-Fock calculations in spherical symmetry. 


approximation Eq. (l34l >. are shown in the case of the SLy4 
interaction in the upper panel of Figure [7] We can see that 
the mass dependence of the diffuseness parameter a in the 
general case is very small. This agrees with the findings of 
Ref. f20J (green dash-dotted lines), where the diffuseness pa¬ 
rameter was extracted from a fit of Hartree-Fock density pro¬ 
files. Considering only the surface term we get a « 0.45 fm, 
while we can observe that taking into account terms beyond 
surface (curvature and mass independent), the diffuseness is 
shifted to lower values of the order of a ss 0.4 fm. This rela¬ 
tively large effect is due to the fact that the non-local curva¬ 
ture term does not contribute to the diffuseness (see Eq. (OTb ). 
Therefore the effect of the curvature energy is to increase the 
local component, which tends to favor a low diffuseness. 

The energy per nucleon is shown in the lower panel of 
Fig. [3 for the three models considered in the upper panel, 
and in comparison to HF calculations. We can see from this 
figure that the variational approach systematically produces 
more binding than the use of a fitted value for the diffuseness, 
as we could have anticipated. Indeed the value of Ref. (20] 
was obtained from a fit of the density, which does not guaran¬ 


tee a minimal energy. Less expected is the fact that the ener¬ 
gies calculated with the three different choices for the diffuse¬ 
ness are very close, though the value of the diffuseness are 
quite different. Specifically, implementing the different dif¬ 
fusenesses into Eq. Q. the resulting total energy reproduces 
the Hartree-Fock nuclear energies with very similar accuracy. 

We can then conclude that introducing higher order terms 
in the variational derivation of the diffuseness, as it has been 
done in equation (l35l >. does not significantly improve the pre¬ 
dictive power of the model. Therefore we will preferentially 
use the simpler expression of the slab diffuseness given by 
equation ( I34| i. This choice is made in all the following fig¬ 
ures, unless explicitly specified. 

C. Decomposition of the surface energy 



FIG. 2. (Color online) Numerical (black circles) and analytical 
(full red line) surface energy per nucleon (see text), and its analytical 
decomposition into plane surface ( °c A 2 / 3 , dashed-dotted blue line), 
curvature (A 1 / 3 , dotted green line), and mass independent (double 
dotted black line) components (eqs j30t . 011 . 03) of symmetric 
nuclei, as a function of the mass number. 

In this section, we study the functional behavior of the an¬ 
alytical formulas of section IIIB I For these applications, we 
keep on focussing on a specific Skyrme interaction, namely 
SLy4 jH. 

In order to verify the accuracy of the analytical expression 
for the surface energy E s , we compare in Figure[2]the sum of 
eqs. (I30l i- (l32l > with the numerical integration of Eq. ( 1 1 2k as 
a function of the nucleus mass. We can see that the analytical 
expressions (full red line) very well reproduce the numerical 
values of E s (black circles). An error smaller than 50 keV per 
nucleon is obtained for the lightest considered nuclei, which 
rapidly vanishes with increasing mass. The deviation for light 
nuclei comes from the approximation in the relation between 
the radius R and the mass A. Indeed, the expansion of the ra¬ 
dius parameter Eq. ® leads to an expansion up to A 1 rs (a jr sat ) 
for E s . The missing terms <x A~ 1 / 3 (a/r sat ) 4 rapidly vanish 
with A, explaining the excellent reproduction of the exact nu¬ 
merical integral. 
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Figure [2] also shows the plane surface, curvature and 
A-independent energy per nucleon components defined in 
eqs. (f30b . (f3Tb and (1321 ). Comparing the total surface energy 
E s (full red line) with E sur f (dashed-dotted blue line), we can 
see that the A 2 / 3 dependence dominates over the whole mass 
table. However, the curvature part (dotted green line), which 
represents the energetic cost of a spherical geometry, cannot 
be neglected even for heavy nuclei, impacting the total en¬ 
ergy of > 300 keV per nucleon for the heaviest nuclei. For 
lighter nuclei (A < 100), the curvature contribution to the to¬ 
tal finite-size effects is of the order of ~ 20%. Though the A- 
independent energy (black dotted line) can be neglected from 
A > 100 for which E mr jj A < 50 keV, it should be taken into 
account for light nuclei if high accuracy is requested. Indeed, 
for A = 40, the A-independent term contributes ~ 5% of the 
total surface energy. 



AAA 


FIG. 3. (Color online) Surface energy per nucleon of symmetric 
nuclei using different choices for the diffusivity parameter a. Panel 
a): variational diffuseness including finite size effects from Eg. ( 135b : 
panel b): variational diffuseness neglecting curvature terms from 
Eq.01; panel c): diffuseness fitted from HF calculations in Ref. na 
a = 0.54 fm. Red lines: total surface energy per nucleon. Blue 
(green) lines: local (non-local) part multiplied by two. 

We now turn to the decomposition of the surface energy 
into a local and a non-local component. It was shown in 
Ref. JH that the local and non-local terms are expected to 
be exactly equal in the case of symmetric matter in a semi- 
infinite slab geometry. This result comes from the fact that 
the one-dimensional Euler-Lagrange variational equation can 
be solved by quadrature 0 As a consequence, it is easy 
to show that if the density profile is the exact solution of the 
Euler-Lagrange variational equation, the first moment of the 
Euler-Lagrange equation implies that the contribution of the 
local term in the surface energy density is at each point of 
space equal to the contribution of the non-local term, lead¬ 
ing to the global equality between the local and non-local slab 
surface tensions: 

a L = a NL . (36) 

Extended to finite nuclei, this result would imply that only 
the local properties of the interaction (that is: the equation 


of state) are needed to predict the surface properties of finite 
nuclei. 

In this paper, we do not solve the Euler-Lagrange equation 
since we impose a given density profile, but we do use a varia¬ 
tional approach in minimising the energy to obtain the diffuse¬ 
ness parameter. Therefore, it is easy to show that our model 
verifies the previous theorem in the one-dimensional case. In¬ 
deed, using the slab diffuseness Eq. (l34k equation (124} reads. 


= a NL = lim - 


1 E 


slab 


°o 2 A 2 / 3 



At first sight this result might look surprising since we have re¬ 
duced the full variational problem to the variation of a single 
variable, which represents a very poor variational approach. 
Equality ( 137} simply means that verifying the Euler-Lagrange 
first moment is equivalent to minimising the energy with re¬ 
spect to a single free parameter. That is, the density derivative 
is well described by the same parameter, here the diffuseness 
a, as the density itself. 

Unfortunately, this elegant theorem cannot be extended to 
the case of a spherical geometry. Indeed, it is easy to show 
that the integrated Euler-Lagrange first moment leads to [25] 


pNL 

°*s 



£NL{r') 


(38) 


The addition of this non-zero integral to the local energy is due 
to the gradient part (°= 1/r) of the spherical Laplacien, which 
comes from the difference between the plane and the spheri¬ 
cal geometry, that is the spatial curvature. Eq. (138} shows that 
in a three-dimensional geometry the equality between the lo¬ 
cal and non-local terms is violated for all components of the 
surface energy, including the term °c A 2 / 3 . 

The left panel of Figure[3]displays the decomposition of the 
surface energy between local (dashed-dotted blue line) and 
non-local (dotted green line) components, when the diffuse¬ 
ness of the density profile is consistently obtained from the 
numerical solution of the variational equation Eq. (135} . We 
can see that the two terms are indeed different. This differ¬ 
ence is however small, and the non-local energy only slightly 
dominates over the local one. This difference is amplified if 
the ansatz for the density profile deviates from the variational 
one. As an example, the central panel in Figure [3] shows the 
surface energy obtained if the simpler expression Eq. (134} for 
the diffuseness is employed. The diffuseness extracted from 
a numerical fit of Hartree-Fock density profiles is employed 
following 0 in the right panel. We can see that the dif¬ 
ference between local and non-local terms is increased as we 
consider density profiles increasingly deviating from the exact 
variational result. 

As we have already remarked, a higher diffusivity (from a) 
to c)) trivially leads to a globally higher surface energy. More 
interesting, the increased deviation from the exact variational 
result from a) to c) leads to a considerable increase of the local 
energy over the non-local one. This is a direct consequence of 
Eqs. @-(l32}. 

From Eq. (|38} . it is clear that the degree of violation of 
equality ( 136} will depend on the functional, as well as on 
the variational model. This point is illustrated in Figure [4] 
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FIG. 4. (Color online) Hartree-Fock calculations. Surface en¬ 
ergy per nucleon (red stars) and its local (blue circles) and non-local 
(green squares) components multiplied by 2, for symmetric nuclei, as 
a function of the mass number. Left (right) panel: Coulomb energy 
excluded (included). 


which shows again the decomposition of the surface energy 
E s into local (blue circles) and non-local parts (green squares), 
calculated numerically from spherical Hartree-Fock calcula¬ 
tions. In the calculations presented in the left panel the 
Coulomb energy, which breaks the equality = E^ L even in 
one-dimensional matter (25tl . is artificially switched off. We 
can see that the Euler-Lagrange result in the slab geometry 
Eq. ( l36l > is reasonably well verified within 10%, especially for 
medium-heavy nuclei A > 90. This shows that the approxi¬ 
mate equality between local and non-local terms is not limited 
to the ETF variational principle, but it is also verified by the 
Hartree-Fock variational solution. However, if the Coulomb 
interaction is included (right panel), the self-consistent modi¬ 
fication of the Hartree-Fock density profile due to Coulomb is 
sufficient to lead to a strong violation of the equality between 
local and non-local terms, going up to 50%. 

This discussion shows that the exact shape of the density 
profile, and in particular the exact value of the diffuseness pa¬ 
rameter, are not important for the determination of the global 
surface energy, but are crucial for a correct separation of local 
and non-local components. In practice it is very difficult to ex¬ 
tract precisely the diffusivity coefficient from theory or experi¬ 
ment: as we have seen in Fig. [I] the diffuseness extracted from 
the Hartree-Fock variational density profile is very different 
from the ETF value, though the energies are close. Moreover 
the equality theorem is violated both because of curvature ef¬ 
fects and of isospin symmetry breaking terms which cannot be 
neglected even in symmetric nuclei because of the Coulomb 
interaction. For all these reasons, we conclude that the con¬ 
tribution from non-local terms cannot be estimated from the 
local part making use of Eq. < [36t . As a consequence, nuclear 
surface properties cannot be understood without mastering the 
gradient and spin-orbit terms of the energy functional. 


HI. ASYMMETRIC NUCLEI 

We now turn to examine the general problem of an ETF 
analytical mass model for asymmetric nuclei, which requires 
the introduction of the proton and neutron density profiles as 
two independent degrees of freedom. In this general case, the 
ETF energy integral cannot be evaluated analytically. The 
usual approach in the literature consists in calculating the 
integral numerically,_with density profiles which are either 
parametrized lf9l JT5L l20ll . or determined with a variational cal¬ 
culation E [26 h 28II. The limitation of such approaches is 
that the decomposition of the total binding energy into its 
different components (isoscalar, isovector, surface, curvature, 
etc.) out of a numerical calculation is not unambiguous nor 
unique 0. Moreover, a numerical calculation makes it hard 
to discriminate the specific influence of the different physi¬ 
cal parameters (EOS properties, finite range, spin orbit, etc) 
on quantities like the surface symmetry energy or the neutron 
skin. 

As a consequence, correlations between observables and 
physical parameters requires a statistical analysis based on a 
large set of very different models. In this way, one may hope 
that the obtained correlation is not spuriously induced by the 
specific form of the effective interaction 112911 . The correla¬ 
tion may also depend on several physical parameters and the 
statistical analysis becomes quite complex [30]. 

Earlier approaches in the literature have introduced ap¬ 
proximations in order to keep an analytical evaluation pos¬ 
sible 0. These approximations however typically neglect 
the presence of a neutron skin, and more generally of inho¬ 
mogeneities in the isospin distribution J9|]. As a consequence, 
the results are simple and transparent, but their validity out of 
the stability valley should be questioned. 

One of the main applications of the present work concerns 
the production of reliable mass tables for an extensive use in 
astrophysical applications 0. For this reason, we aim at 
expressions which stay valid approaching the driplines. In 
the specific application to the neutron star inner crust, even 
more exotic nuclei far beyond the driplines are known to be 
populated (33 . H.3;]. We will not consider this situation in 
the present paper, because a correct treatment of nuclei be¬ 
yond the dripline imposes considering the presence of both 
bound and unbound states which modify the density pro¬ 
files and leads to the emergence of a nucleon gas. Opti¬ 
mal parametrized density profiles have been proposed for this 
problem fl3l20l . [34}l . but the developement of systematic ap¬ 
proximations to analytically integrate the ETF functional in 
the presence of a gas is a delicate issue, which will be ad¬ 
dressed in a forthcoming paper |25i]. 

A. Decomposition of the nuclear energy 

The presence of two separate good particle quantum num¬ 
bers, N and Z, implies that we have to work with a 2- 
dimensional problem, and introduce, in addition to the total 
density profile Eq. ©, an additional degree of freedom. Con¬ 
cerning the energy functional, it is customary to split it into an 
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isoscalar and an isovector component: 

'#’{p,p3\ = [p, p3 = o] + JO™ [p, p 3 ] (39) 

with: 

ft 2 

Jf IS [p,P?,] = —T + C eff pT+{C 0 + C 3 p a )p 2 

+ Cfi,i (Vp)" + C so J • Vp, (40) 

N [P , P3] = ^ [P, P3] - ^ [P, P3 = 0] 

+ D effP3^3 + {Do + D 3 p a )p 2 
+ Df in (Vp 3 ) 2 + D so i 3 ■ Vp 3 , (41) 

where we have introduced the local isoscalar and isovector 
particle densities, kinetic densities and spin-orbit density vec¬ 
tors. Isoscalar densities are given by the sum of the corre¬ 
sponding neutron and proton densities, while isovector den¬ 
sities (noted with the subscript ”3”) are given by their differ¬ 
ence. As for symmetric matter, the semi-classical Wigner- 
Kirkwood development in ft allows expressing all these den¬ 
sities in terms of the local isoscalar p = p n + p p and isovector 
P3 = Pn ~ Pp density profiles, as well as their gradients. In 
equation d40l >. the isoscalar energy density also depends on p 3 
because of the presence of the kinetic densities x = X n + x p 
which cannot be written as a function of p only. Therefore, to 
truly obtain the isoscalar part in Eq. (l39t , we have to consider 
J4? IS [p,p 3 = 0]. The iso-vector energy density Eq. ( 14 1 t i con¬ 
tains therefore terms which explicitly depend on the isovector 
densities, but also an isovector contribution of the so-called 
isoscalar component ,W, AS . Detailed expressions, and defini¬ 
tion of parameters are given in appendix [A] 


1. Isospin inhomogeneities 


Concerning the density profiles, we choose to work with the 
total density p(r) and with the proton density profile p p (r). 
Alternatively, we could as well have used (p,p„) or (p p ,p„) 
as independent variables, and we have checked that these 
different representations lead to the same level of reproduc¬ 
tion of full Hartree-Fock calculations. The total density is 
parametrized by Eq. ©. where now the saturation density pa¬ 
rameter p sat corresponds to the equilibrium density reached 
in asymmetric matter 0. This density depends on the asym¬ 
metry 5 which represents the nucleus bulk asymmetry, defined 
below: 


Psat (8) = psat { 0) 


3 L sym 8 2 \ 

K$at 4 ” K S y m 8 2 J 


(42) 


In this expression, K sat = 9p 2 at d 2 (Jf? / p)/dp 2 \ Psal 
is the nuclear (symmetric) matter incompressibil¬ 
ity, and Lsym — 8 psat c) (■^sym j P)f 8 p \ p sfl! and K S y m = 
9p 2 a td~(^sym/P)/dp~ | Psat are the slope and curvature of 
the symmetry energy at (symmetric) saturation, where we 
have introduced the usual definition of the symmetry energy 
density : 


'ffisym — ~ P 


1 i d 2 J$? 


dp 2 


P3= 0 


(43) 


As a consequence, the radius parameter/? entering Eq. © also 
depends on the nucleus bulk asymmetry 5. Indeed, in Eq. ®, 
the equivalent homogeneous sphere radius now reads Rhs = 
A i / 2 ’r sat {8), where the mean radius per nucleon is r sat (8) = 

(lnp sat (8)y 1/3 . 

The proton density profile is parametrized as an indepen¬ 
dent Fermi function [20): 

P P (r) = P sat, P F p (r) ; F p (r) = (l +e (^)/«p) _1 . (44) 

In equation (l44t . the proton radius parameter R p is deter¬ 
mined, similarly to Eq. ©, by proton number conservation 
as: 


R p = R/isp 





,(45) 


with Rhs p (8 ) = / 2 r sa t p(8) the equivalent homogeneous 

proton sphere radius, r sa t :P (8) = (%np sa t, p (8)) and 
where we assumed a p ^R p . 

The diffusenesses a and a p will be calculated in sec¬ 
tion IIIIBI by a minimization of the surface energy, as it has 
been done for symmetric nuclei in section HTbI where a p = a. 
We can anticipate that the isoscalar diffuseness a will be mod¬ 
ified with respect to the result of symmetric nuclei Eqs. (l34ll 
and (l35l >. 

In order to have the correct bulk limit of infinite asym¬ 
metric matter, the parameters p sat and p m i. P introduced in 
Eqs. © and (l44l) respectively represent the saturation densi¬ 
ties of baryon and proton of asymmetric matter. These densi¬ 
ties are related to the properties of the Skyrme functional and 
to the bulk asymmetry 5 = 1— 'Ipsat , P IPsat by Eq. (l42l) . 

The bulk asymmetry differs from the global asymmetry 
I = 1 — 2Z/A because of the competing effect of the Coulomb 
interaction and symmetry energy, which act in opposite di¬ 
rections in determining the difference between the proton and 
neutron radii j26ll27l [35il : 


8 = 


i i 3«c Z 2 
1 + 82 J573 


1 + 


9 J„- 


1 


4Q Jpi 


(46) 


In this equation, J sym = <7t?sym[Psat\/Psat is the symmetry en¬ 
ergy per nucleon at the saturation density of symmetric matter, 
Q is the surface stiffness coefficient, and a, is the Coulomb pa¬ 
rameter. Because of the complex interplay between Coulomb 
and skin effects, the bulk asymmetry 8 of a globally sym¬ 
metric 7 = 0 nucleus is not zero, though small for nuclei in 
the nuclear chart. We have shown in Ref. El that account¬ 
ing for the 8 dependent saturation density gives a reasonably 
good approximation of the isospin symmetry breaking effects 
in I = 0 nucleL_A complete discussion on this point can be 
found in Ref. lf36h . 

As a consequence, the interval of 8 is slightly smaller than 
the interval of I over the periodic table. The relation between 
the global asymmetry and the asymmetry in the nuclear bulk is 
shown in Fig.© From this figure we can see that 5 is a slowly 
increasing function of the global asymmetry I. This value in¬ 
creases to —0.1 < 8 < 0.3 if we consider the ensemble of the 
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2. Bulk energy: limit of asymmetric nuclear matter 


Following the same procedure as for the symmetric case, 
we can define the bulk energy in asymmetric matter as: 

E b {8) = J%at{8)V HS {8) = A rar (5)A, (47) 
J%a,(8) — ^sat (■ 8)Psat(8 ), (48) 

where Vhs(8) = 4/3 7iR 3 HS (8) =A/p sat (8) is the equivalent 
homogeneous sphere volume and A sa t(8) corresponds to the 
chemical potential of asymmetric nuclear matter: 




dp 


lPsat(8),p sali3 (8)] 


= A sat(8) 


P 


■ (49) 

[Psal(8)iP sa t,3(8)} 


FIG. 5. (Color online) Bulk asymmetry Eq. <42} as a function of the 
global asymmetry I for nuclei within the theoretical driplines eval¬ 
uated from the SLy4 energy functional. The different colors corre¬ 
spond to different intervals in mass number: 40 < A < 100 in red, 
100 < A < 150 in blue, 150 < A < 200 in green, A > 200 in grey. 
The function y = x is also plotted (black). 


heavy and medium-heavy nuclei within the driplines [!;65]]. It 
is also observed from Fig. [5] that as the mass A increases, 5 
becomes closer to /, as expected from the analytical expres¬ 
sion <46}. 



N 

FIG. 6. (Color online) jS-stable nuclei (green), unstable nuclei syn- 
thetized in the laboratory 01 (red), theoretical neutron and proton 
driplines evaluated from the SLy4 energy functional (black squares) 
and some iso-5 lines (blue dots) are plotted in the N,Z plane. 


Figure [6] shows in the ( N,Z ) plane the heavy and medium- 
heavy measured nuclei, the theoretical neutron and proton 
driplines evaluated from the SLy4 energy functional, and 
some iso-5 lines. We can see that all A < 40-isotopes ever syn¬ 
thesized in the laboratory lay between 5 « 0 and 5 ps 0.2. Fur¬ 
thermore, the theoretical neutron dripline well matches with 
the iso-5 line 5 ps 0.3, which roughly corresponds to / ft; 0.4. 

This means that in the following, we will be interested in 
approximations producing reliable formulae up to 5 ps 0.3. 


Here, p sa tp = Psat ~ 2-Psat,p, and the total energy density Jff is 
given by Eq. <39}. 


3. Decomposition of the surface energy 


The surface energy E s (8) corresponds to finite-size effects 
and can be decomposed, as in the symmetric case in section HI1 
into the plane surface, the curvature, and the higher order 
terms. It is defined as the difference between the total and 
the bulk Eb(S) energy. 


E s (8) = J drJflp,p 3 ]-Ji? sat (S)V HS (8) 

poo 

= 4 nj dr{jf[p,p 3 ]-A^(5)p}r 2 . (50) 


Because of the isospin asymmetry, the Skyrme functional 
now depends on the two densities p and p 3 = p — 2 p p and on 
the two gradients Vp and Vp 3 = Vp — 2Vp p . 

Making again the decomposition of the energy density into 
an isoscalar (only depending on the total density) and an 
isovector component (depending on p and p 3 ), we get from 
Eqs. <40} and (|4TT) : 



E s = E\ 

with 


II 

E? 

hf 

^ dr|^ /5 [p,p 3 = 

=*’J 

^ dr|^T[p,p 3 = 

•*t 

II 

hf 

^ d r^Jt IV \p,p 3 ] 

= 4n ) 

f drjjT[p,p 3 ]- 


•is , E iv, 


(51) 


nj _ ^ IS \P.satiPsat, 3 — 0 ] 1 2 

‘ J Psat(S) P J 

nj [psat, Psat,3 = 0] 

J PUS) 

^ IV \p sutl p sa0 ] 

Psat (8) 

■ : yP \p sa t ■ Psa! 


p 

(52) 

2 


Psat (8) 


P >r 


P 


(53) 


It is interesting to remark that Eq. <50} is not the only pos¬ 
sible definition of the surface energy in a multi-component 
system. Indeed in a two-component system, there are two 
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possible definitions of the surface energy which depend on 
the definition of the bulk energy in the cluster 11381 - 44011 : the 
first one is given by Eq. (f50b and corresponds to identify¬ 
ing the bulk energy of a system of N neutrons and Z pro¬ 
tons to the energy of an equivalent piece of nuclear matter . 
The second definition E S =E — jl n N — fi p Z + pV corresponds 
to the grandcanonical thermodynamical Gibbs definition, and 
gives the quantity to be minimized in the variational calcu¬ 
lation conserving proton and neutron number. Though this 
second definition has been often employed in the ETF liter¬ 
ature [HI IDj - M . the first one Eq. ( l50l > is the most natural 
definition in the present context. Indeed, using the decomposi¬ 
tion Eq. (l39l > between isoscalar and isovector energy densities, 
only this definition allows recovering for the isoscalar energy, 
the results of section [IT] concerning symmetric matter. More¬ 
over, we have shown in Ref. IHill that the best reproduction of 
full Hartree-Fock calculations is achieved considering that the 
bulk energy in a finite nucleus scales with the bulk asymmetry 
8 as in Eq. (l50l >. rather than with the total asymmetry I, as it 
is implied by the Gibbs definition. 

Let us first concentrate on the isoscalar surface energy. The 
dependence of the surface energy on the bulk asymmetry 5 
implies that its decomposition into an isoscalar and an isovec¬ 
tor part is not straightforward. Indeed, although the isoscalar 
energy E 1 / does not depend on the isospin asymmetry profile 
P 3 (r), it does depend on the bulk isospin asymmetry 8 through 
the isospin dependence of the saturation density p sa t(8) ap¬ 
pearing in the density profile p Eq. ©■ Moreover, in Eq. < l32b 
the isoscalar bulk term which is removed depends directly on 
8 too, because of the equivalent volume Vhs = A / p sat (8). The 
quantity E 1 / has therefore an implicit dependence on isospin 
asymmetry 5. 

The isoscalar surface energy can be calculated exactly for 
any nucleus of any asymmetry, with the expressions devel¬ 
oped in section QI] In particular we can distinguish a plane 
surface, a curvature, and a mass independent term: 

4 s =KU+ E, i n +^+o ((^§) <A ~ 1/3 ) < 54 > 

with an identical result as in Eqs. (l30l >. (OH . ED . namely: 


pis 

^ surf 


E IS = 

^curv 


t?IS _ 

t ind ~ 


c/?L 
®surf ' 


°curv 




1 


C/3‘ 


>NL 


a 2 (A,S) surf 

1 


a 2 {A,8) 


1 


a 2 (A, 8) 


aeNL 

®curv 


^NL 

°ind 


a(A,8) 
r S at{8) 
a(A,8) 
r sa t{8) 
a(A, 8) 
rsat(S) 


a 2 /\ 


A 1 / 3 , 


(55) 

(56) 

(57) 


The local and non-local functions are given by 
Eqs. ( fl9b and (l22l >. where the saturation density now depends 
on asymmetry p sat = p sat (8 ) through Eq. (l42l >. The other dif¬ 
ference with respect to the case of symmetric nuclei Eqs. ( l30l ). 
(ED, (ED. is that now the diffuseness depends on the asym¬ 
metry 5. 

Since the analytical expressions of the isoscalar surface en¬ 
ergy E[ s are the same as in symmetric nuclei, the same accu¬ 
racy and conclusions as in section QI| are dressed: we can vari- 



FIG. 7. (Color online) Diffuseness as a function of the isospin 
asymmetry, for four isobaric chains (A = 400: full lines, A = 200: 
dotted lines, A = 100: dashed-dotted lines, A = 50: dashed lines). 
Red lines: calculations using the slab diffuseness Eq. l l34l i. Blue 
lines: calculations using the spherical diffuseness Eq. ED Green 
lines: calculations using the quadratic diffuseness, fitted from HF 
density profiles in Ref. [2(1]. 


ationally evaluate the isoscalar diffuseness a, solving equa¬ 
tion ED. or using equation ED which amounts to neglect¬ 
ing terms varying slower than A 2 / 3 . Though we have con¬ 
sidered only isoscalar terms, the diffuseness a does depend 
on the isospin asymmetry 8 because of the 5 dependence 
of the saturation density. These results, as well as the fit 
from HF density profiles [20!], where mass independence and 
quadratic behaviour in 8 is assumed (that is: a = C\ + CAS 1 ), 
are shown in Fig. [7] Concerning the mass-dependence of 
Eq. ED (blue lines labelled ”Eq. ED ), we observe a slight 
spread for masses from A = 50 to A = 400, corroborating both 
the mass independence assumption in the HF fit j20il and the 
previous conclusions in section IU bI to obtain the diffuseness 
we can neglect the mass dependence and limit to terms A 2 / 3 
(red line, labelled ”Eq. ED”)- However, one can see that the 
dependence found from the variational equation is opposite 
to the one exhibited by the fit to HF results: the diffuseness 
decreases with 5 instead of increasing. It is difficult to be¬ 
lieve that such a huge and qualitative difference might come 
from the difference between ETF and HF. The discrepancy 
rather suggests that the variational procedure should include 
the isovector energy to obtain the correct behaviour of the dif¬ 
fuseness with the isospin asymmetry. Indeed, we will see in 
section lHIBI that adding the isovector part reverses the trend. 

This discussion shows that, in the case of asymmetric nu¬ 
clei, Eq. ED which only takes into account the isoscalar 
terms, is not a good approximation to find the diffuseness. 
This statement is confirmed by Fig. [8] where the isoscalar sur¬ 
face energy per nucleon is plotted for different isobaric chains 
and for different prescriptions for the diffuseness. The full 
red and the dashed-dotted lines stand for the diffuseness given 
by Eq. ED and Eq. ED respectively. There is almost no 
difference in the isoscalar surface energy for these two pre¬ 
scriptions. In addition, the observed 8 dependence is ex- 
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FIG. 8. (Color online) Isoscalar surface energy per nucleon as a 
function of the isospin asymmetry, for four isobaric chains. Full red 
lines: calculations using the slab diffuseness Eq. tU- Dash-dotted 
blue lines: calculations using the spherical diffuseness Eq. d35). 
Dashed green lines: calculations using the quadratic diffuseness, fit¬ 
ted from HF density profiles in Ref. [20t l. 


tremely weak. The isoscalar surface energy evaluated with 
the quadratic diffuseness [[ 20 ] is represented in dashed green 
line. A qualitative and quantitative difference is observed with 
respect to the two other curves. This indicates again that the 
isoscalar and isovector component of the surface energy can¬ 
not be treated separately, and the correct 8 dependence of the 
isoscalar surface energy, as well as of the isoscalar diffuse¬ 
ness, requires to consider the total surface energy in the varia¬ 
tional principle. 

It is also interesting to analyse the 8 dependence of the sur¬ 
face symmetry energy based on the fitted quadratic diffuse¬ 
ness: its sign is positive, which contrasts with studies based 
on liquid-drop parametrizations of the nuclear mass HE 
I35[l4ll - l43ll . This behavior is due to our choice of definition of 
the surface in a two component system, as discussed at length 
in Ref. lfl9l] . 


B. Approximations for the isovector energy 

In this section, we focus on the residual isovector surface 
part defined by Eq. (|53| >. which cannot be written as inte¬ 
grals of Fermi functions as in the previous sections. Indeed, 
the isovector density P 3 appearing in the energy density is not 
a Fermi function, meaning that it cannot be analytically in¬ 
tegrated to evaluate E™. Approximations are needed to de¬ 
velop an analytic expression for this part of the energy, and 
we will consider in the following two different approaches. 
At the end, we will verify the accuracy of our final formulae, 
in comparing the analytical expressions with HF calculations. 


1. No skin approximation 

As a first approximation, we neglect all inhomogeneities 
in the isospin distribution in the same spirit as Ref. 0 . 
This simplification consists in replacing the isospin asym¬ 
metry profile pT,(r)/p(r) in Eq. (l53l > by its mean value (8). 
Within this approximation, the local isovector energy only de¬ 
pends on the total baryonic density profile p defined Eq. <[ 6 ]), 
and the non-local isovector part, involving gradients Vp 3 , is 
identically zero. In other words, this approximation amounts 
neglecting the non-local contribution to the isovector surface 
energy. 

Integrating in space the equality Pi(r) = (S)p(r) we imme¬ 
diately obtain that the mean value of the isospin distribution 
is given by the global asymmetry of the nucleus: 

(8) = — = I- ( 58 ) 

In particular, in this approximation, the bulk isospin asym¬ 
metry 8 is equal to the global asymmetry /, at variance 
with the more elaborated relation between 8 and I given by 
Eq. d46l >. In neglecting isospin inhomogeneities, we indeed 
neglect both neutron skin and Coulomb effects which are re¬ 
sponsible for the difference between 8 and I. Consequently in 
this section, the saturation density p sat of asymmetric matter is 
still given by Eq. d42l) . but replacing 8 by I. This no-skin ap¬ 
proximation therefore modifies the bulk energy Eq. (l47l) . and 
the isoscalar energy Eq. (fl2l ). 

The choice of I instead of 5 to compute the saturation den¬ 
sity only slightly worsens the predictive power of the total 
ETF energy with respect to Hartree-Fock calculations, but the 
relative weight between bulk and surface energies is drasti¬ 
cally modified. In particular, this change of variable switches 
the sign of the symmetry surface energy 0 ]. 

The obvious advantage is that analytical results can be ob¬ 
tained without further approximations than the ones devel¬ 
oped in section iHAl as we now detail. 

Replacing p 3 (r) by Ip{r) and p sat p(8) by Ip sa ,(I) in 
Eq. (l53b . allows to express the energy density as a function 
of p (r) only. Thus we can follow the same procedure as for 
symmetric nuclei in section UTAl and analytically integrate the 
energy density. Making a quadratic expansion in I for the ki¬ 
netic densities T 3 gives the following expressions: 



where X^ y = {C e //, a,D e ff } stands for the effective interac¬ 
tion parameters appearing in the isovector local terms. The 
coefficients i ff v are given by: 
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(62) 


where m/m* at = {m/m* at n + m/in* sat p )/ 2, 8m sat = ( 8m sat ,n + 
Sm S at,p)/ 2, Am raf)3 = (m/m* atn - m/m* atp )/{21) = 

{8m sa t,n — 8m sat p)/(21), and where the coefficients Py k> 
are defined by equation ( IB2b . 

As for the isoscalar energy, Eq. d59l ) shows that the dom¬ 
inant finite-size effect is a surface term (°= A 2 / 3 ). Additional 
finite-size terms, which would be absent in a slab configu¬ 
ration, are found in spherical nuclei. As we have only con¬ 
sidered the local part of the isovector energy, we recover the 
same diffuseness dependence as in the local isoscalar terms 
Eqs. (fl5l i and (IT6l) . 

We have seen in section IHBI that the diffuseness a can be 
obtained by minimization of the energy per nucleon with re¬ 
spect to its free parameters. In this no-skin approximation, the 
only non-constrained parameter of the model is again the dif¬ 
fuseness parameter a, as for symmetric nuclei. Therefore, we 
can apply Eq. ( l33l ) in order to obtain the ground state energy. 
If we neglect the curvature and mass independent terms, we 
obtain an expression similar to Eq. ( l34t : 


a = 



\ 


KrfW + KrfW 21 


(63) 



I 

FIG. 9. (Color online) Diffuseness as a function of the global asym¬ 
metry, for four isobaric chains (A = 400: full lines, A = 200: dotted 
lines, A = 100: dashed-dotted lines, A = 50: dashed lines). Red 
lines: calculations using the slab diffuseness Eq. <631 . Blue lines: 
calculations using the spherical diffuseness Eq. i[64]i. Green line: 
calculations using the quadratic diffuseness fitted from HF density 
profiles 12011. Grey line: calculations using the diffuseness Eq. <661 . 
based on 1161. 


where the coefficients ^‘urf depend on the saturation density 
Psat(I)- This expression corresponds to the diffuseness of one¬ 
dimensional semi-infinite asymmetric matter. Considering all 
the terms of Eq. (f59l> . the diffuseness corresponding to the 
complete variational problem is given by the solution of the 
following equation: 

3 ^ tul /“ ) ( -+ 2 furv + ^ mrv^ ) A1 ^ 1 

+ ('K 1 „/ + *2P ! )A 2/3 + 2_*;»A (JlY (64) 

\ r sat / \ r sat J 

-^C/A 2/3 = 0. 

r sat 

Figure[9]displays the results of Eqs. ( |63| > and ( l64t . At vari¬ 
ance with Fig.[7]where we only took into account the isoscalar 
energy, we can clearly see that adding the isovector energy to 
the variational procedure leads to the expected behavior of a 
diffuseness increasing with asymmetry. 

This behavior shows the importance of the isovector part to 
correctly determine the isoscalar diffuseness a. As for sym¬ 


metric nuclei, we observe again that the mass dependence of 
the diffuseness calculated in the spherical case is negligible 
(only a slight spread of the blue curves, no spread in the red 
curves). 

The analytical total surface energy E s = E 1 / + E[ v per nu¬ 
cleon, given by Eqs. (IT5l) . (IT6l) and (l59l >. is plotted on Fig.flOl 
for different isobaric chains. The results using the slab dif¬ 
fuseness (full red curves) are very close to the ones obtained 
by solving Eq. (l64l) (dash-dotted blue curves), and to the ones 
using the numerical fit to HF calculations of Ref. 0 (dashed 
green curves), even if the corresponding values for the a pa¬ 
rameter are very different. The conclusions are thus the same 
as in section IHBI although curvature (and mass independent) 
terms are important to reproduce the energetics, they are not 
required to determine the diffuseness. Therefore this latter can 
be well determined by the simplest expression, Eq. (l63l) . 

For completeness, we also compare our results to the ap¬ 
proximation for the surface energy proposed in Ref. [16], and 
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FIG. 10. (Color online) Total surface energy per nucleon as a func¬ 
tion of the global asymmetry, for four isobaric chains. Full red lines: 
calculations using the slab diffuseness Eq. ( I63l >. Dash-dotted blue 
lines: calculations using the spherical diffuseness Eq. m- Dashed 
green lines: calculations using the quadratic diffuseness fitted from 
HF density profiles li201. Dotted grey lines: calculations using the 
diffuseness Eq. ( 1661 . based on 01. 


represented by grey curves in Figs.l9landlT0l 


E s =£'/(! = 0) 


£f (/ = 0) L sym 
AV 3 K sat 


a 




fsat (J — 0 ) 


A 2 / 3 I 2 .( 65) 


In Ref. 10 ], no expression for the diffuseness was proposed. 
For consistency, we have determined the a parameter entering 
Eq. ( |65| > by minimizing the surface energy given by the same 
equation, leading to: 


Krf{I = 0) 

I 1 i 2L svm t 2 ) 

i, 1 + K sal 1 J 

_ f 

^W( /= °) 


— 2 

Lsym 

Ksym 

-ttJ 

i 2 


To obtain Eq. (l65l > , the authors of Ref. 0 did the same ap¬ 
proximation P 3 (r) = Ip (r) as we made, neglected the curva¬ 
ture and constant terms, and assumed the equality E[: = E^ L 
for the isovector part in order to evaluate the non-local isovec¬ 
tor energy. As we have shown in section III Cl this property 
fails in a three-dimensional system. As a consequence, the 
diffuseness which is determined by the balance between local 
and non-local parts, is overestimated (see Fig. [9| and finally 
leads to a largely underestimated energy, as seen in Fig.lTOl 
To check the accuracy of our analytical no-skin expression 
given by Eqs. (O, (0) and <[63]», we will quantitatively com¬ 
pare our analytical results with Hartree-Fock calculations in 
section lHIB 31 


2. Gaussian approximation 

To take into account isospin inhomogeneities, we develop 
in this section an alternative gaussian approximation to the 


isovector surface energy. In particular, as in section HIlAI we 
will distinguish the bulk asymmetry 8 Eq. (l46l) from the global 
one I, which allows considering skin and Coulomb effects. 
This approximation is therefore expected to be more realistic 
than the no-skin procedure developed in section UlIB II 



r (fm) 


r (fm) 


FIG. 11. (Color online) Numerical isovector energy density profile 
(red full lines) and Gaussian approximation Eq. ( 168b (black dashed- 
dotted lines) for two masses A = 50 (left curves of each panel) and 
A = 200 (right curves of each panel), a) 8 = 0.1; b) 8 = 0.2; c) 
8 = 0.3; d) 8 = 0.4. 


Since E[ v is the surface isovector energy, the corresponding 
energy density 

■xfM = ^' r M - (67) 

is negligible at the nucleus center, where p -A p sat . This is 
shown in Fig. HU which displays this quantity for several nu¬ 
clei in a representative calculation using the diffusenesses a 
and fl p from Ref. [200, and with the interaction SLy4. More¬ 
over, as it is a surface energy, the maximum is expected to be 
close to the surface radius R, that is the inflection point where 
p(R) = Psat (5)/2. Thus we approximate the isovector energy 
density by a Gaussian peaked at r = R: 

Jf?/ v (r)~& tot (r) = £/(A 1 8)ex p ^(A^S) ) ’ (68) 

where &/ is the maximum amplitude of the Gaussian and fr¬ 
its variance at R: 

^/(A,8) = J(f s IV [p(R),p 3 (R)], (69) 

a 2 (A,8) = ~^(A,8)^-^-j . (70) 

Fig. 1 1 llshows the quality of this Gaussian approximation on 
the energy density profile for several nuclei. Each panel corre¬ 
sponds to a different representative value of 8: 8 = 0.1 (upper 
left) corresponds to most stable nuclei (see Fig. [ 6 ]»; medium- 
heavy neutron rich nuclei synthesized in modern radioactive 
ion facilities lay around 5 = 0.2 (upper right); the (largely 
unexplored) neutron drip-line closely corresponds to 8 = 0.3 
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(lower left); the higher value 8 = 0.4 (lower right) is only 
obtained beyond the dripline, that is for nuclei which are in 
equilibrium with a neutron gas in the inner crust of neutron 
stars. 

We can see that for all these very different asymmetries, the 
exact energy density (full red lines) is indeed peaked at the 
equivalent hard sphere radius R. However, we can notice that 
the profiles have small negative components. We thus expect 
the Gaussian approximation will overestimate the isovector 
energy part. 

As Gaussian functions and their moments are analytically 
integrable, this approximation allows obtaining an analytical 

expression for the isovector energy E[ v tv 4n J r 2 %:,t(r)dr. 

Indeed, neglecting the terms ~ e - fi2 /( 2tT ~) ; we obtain (see ap- 
nendix lUT} : 


e! v 


2(2n) 3 ^ 2 o(A,8,a,a p )g/(A,8,a,a p )r 2 at (8) 
, 2/3 tf 2 (A,<5,a,a p ) 2 n 2 /a(A,<5)\ 2 

r^at(S) 3 \ rsat (S)J 


,(71) 


we need to evaluate the symmetry energy at two different 
densities: at p sa t{8) and at the surface radius where p (R) = 
Psai{8) /2. For this reason, we will apply Eq. ({71} to two dif¬ 
ferent densities p* = p sat ( 0) and p* = p sat ( 0)/2. At p* = 
p sa t( 0), the coefficients (/* . L t . AT* ) are the usual symmetry 
energy coefficients (J S ym,L S y m ^K sym ). Their values for the 
Skyrme interaction SLy4 are J sym = 32 MeV, L sym = 46 MeV, 
and K sym = —119.8 MeV. At one half of the saturation of 
symmetric nuclear matter, p* = p sat {0)/2 we label the corre¬ 
sponding coefficients {J\/ 2 -.L\/ 2) K\ / 2 ) which, for the Skyrme 
interaction SLy4, are J ] / 2 = 22.13 MeV, Li/ 2 = 38.6 MeV, 
and K ] /2 = -74 MeV. 

Using the expansion around p . = p Mf (0)/2 for the first term 
of Eq. (172} and around p* = p sa t{ 0) for the second one, we 
obtain, at second order in 8: 


g/ (A, 8) 
Psat (0) 


J \/2 { AR(a) \ 2 
8 Va(A,5)7 
J 1/2 A R(a) 1 / A R(a) 

2 a(A, 8) 2 \a(A, 8) 


2 


8 


where we have highlighted the dependence on the nuclear 
mass number A, bulk asymmetry 5, and diffusenesses a(A,S) 
and ci p (A,8) when they explicitly appear. The neglected terms 
are of the order (a/r sar ) 4 A~ 2 / 3 . We can notice that the curva¬ 
ture term (°= A 1 / 3 ) is missing. This is due to our approxima¬ 
tion. Indeed, we have assumed that the isovector energy is 
symmetric with respect to the inflection point for which the 
curvature is zero, such that the curvature is disregarded by 
construction. 

Though equation ( 171} is an analytical expression, the ex¬ 
plicit derivation of the amplitude g/(A,8) and of the variance 
<7(A, 8) leads to formulae far from being transparent. In par¬ 
ticular, it is not clear how the different physical ingredients of 
the energy functional (compressibility, effective mass, sym¬ 
metry energy) and of the nucleus properties (neutron skin, dif¬ 
fuseness) affect the isovector surface properties. For this rea¬ 
son, we turn to develop a further approximation for the isovec¬ 
tor energy part E ,v in terms of the nuclear matter coefficients 
J, L and K, and of the neutron skin thickness. Moreover, these 
approximations will allow to find a simple analytical expres¬ 
sion for the diffuseness. 

Making the usual quadratic assumption for the symme¬ 
try energy J^ IV [p, p 2 \ = ^ym[p](p 3 /p) 2 , the amplitude 
gZ(A, 8) Eq. ( 169} reads 

*(A,8) = ^ sym [p(R)} (f^) 2 -^ym[Psat(8)]-A -^§ 2 . 

(72) 


In order to have a simpler explicit expression, we make a den¬ 
sity expansion of the symmetry energy per nucleon e sym [p] = 
M' S ym{p]/p around a density p*, such that: 


't^sym [p] — P 


J - + wZ p - p - )+ m ip - p - f 


,(73) 


where J * — <^sym[p*\/P*> L* — 3p*d(jffsym/p*')/dp\p* and 
K * = 9p 2 d 2 (j4? sym /p*)/dp 2 \ p *. As we can see in Eq. (l72l) . 


/j Jsym \ A R{ci) 

V J i/ 2/ a(A,8) 

1 / L sym L l/2 \ / A R{a)\ 2 
4 { J l/2 K sat ) {a(A,S)J 

Notice that the K sym parameter does not appear in this equation 
because of the truncation at second order in 8. In Eq.(IT4}. the 
isospin asymmetry inhomogeneities clearly appear through 
the quantity A R(a) = R(a) —R p (a) which represents the neu¬ 
tron skin thickness: 

( 7Z 2 Q 2 \ 

AR{a)=/SR HS (\ + ——- - , (75) 

\ d KhsKhs,p J 

where A R HS (A,Z) = AR(a = 0,A,Z) = R HS (A) - R HS , P (Z) is 
the neutron skin thickness of nuclei theoretically described by 
hard spheres. Moreover, we have considered the diffuseness 
difference a — a p as a second order correction with respect 
to the neutron skin, and have assumed a = a p in Eq. (174} . We 
have also used the following expansion in A R(a)/a to evaluate 

P3(*): 

2 Pp(R) = Psa,,p(8) [1 - AR(a)/(2a)] (76) 

+ 0((AR(a)/a) 3 ). 

Eq. ({74} gives a relatively simple and transparent expression of 
the isovector energy density at the nuclear surface, as a func¬ 
tion of the EoS parameters. The situation is more complicated 
for the variance <7(A, 5) which also enters the isovector energy 
Eq.CZD- This quantity involves the second spatial derivative 
of the energy density Eq. (ITO} . therefore its explicit expres¬ 
sion is not transparent, even with the previous simplifications. 
Extra approximations are in order. 

From Fig.HU we can observe that the width of the numer¬ 
ical gaussians, that is the values of cr 2 (A,5), is almost inde¬ 
pendent of the bulk isospin 8. This numerical evidence can 
be understood from the fact that the width gives a measure of 
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the nucleus surface, which is mostly determined by isoscalar 
properties. It is therefore not surprising that the dominant 
isospin dependence is given by the amplitude ,</? which rep¬ 
resents the isovector energy density at the surface. For this 
reason, we evaluate the variance at 5 = 0: 


a(A,8) ss a (A) 


N 


2 


K l/1 

18-6/2 


a o = Oo- 


(77) 


In this equation, stands for the diffuseness at 8 = 0. We re¬ 
call that this quantity does not depend on the nucleus mass if 


7 iv 
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K l/2 

18 ^ 1/2 


r S at(8) 


2 



In principle the surface coefficients {J\/ 2 ,L\/ 2 ,K 1 / 2 ) can be 
expressed as a function of the bulk ones (J sym -, L sym , K sym ) by 
using polynomial expansion in the density. However, we can 
see from Eq. ( I79[ i that the surface isovector energy E[ v is pro¬ 
portional to the symmetry energy J l / 2 evaluated at the sur¬ 
face R. It is quite natural that the surface energy component 
is mainly determined by the surface properties of the nuclei, 
and therefore, the surface symmetry energy is mainly propor¬ 
tional to the isovector parameter J\ / 2 . For this reason, express¬ 
ing Eq. ( 1791 ) only in terms of bulk quantities (7 ivm , L JVm , A " sym ) 
would make Eq. ( f79t less transparent. 

For completely symmetric nuclei, that is A R = 0 and 5=0, 
the isovector energy is identically zero as it should. However, 
if we neglect the neutron skin thickness only, that is we con¬ 
sider A R = 0 but 5 ^ 0, a non-zero isovector surface energy is 
obtained, given by 


j^IV.AR—0 

E surf 



ft Psatify ao 
1 8^1/2 


{Jl/2-Jsym) 5 2 A 2 / 3 . 


(80) 


This expression is proportional to the energy density differ¬ 
ence between bulk and surface (7; n — J sy m) , that is to the 
L sym parameter. In this approximation, the diffuseness a (A, 8) 
does not appear, which means that the isovector surface en¬ 
ergy contributes to the determination of the diffuseness only 
if we consider the neutron skin. 

From a mathematical point of view we can also consider 


we do not take into account terms beyond surface in the vari¬ 
ational approach discussed in section IIIB I This approximate 
mass independence of the variance can be verified in Fig. |TT] 
the width of the two gaussians corresponding to A = 50 and 
A = 200 are very close. Neglecting the isovector component 
at 5 = 0, the diffuseness is then given by the expression ( l34l ) 
valid for symmetric matter: 

ao = y/^f(S = 0)/^ urf (8=0). (78) 

Inserting Eqs. ( 1741 ) and ( 1771 ) into CD, the surface isovector 
energy can be expressed as a function of the symmetry energy 
coefficients (J sym ,L sym 

5 Ksym ) • 


Jsym \ AR[ci) 1 
Xfi)~ a(A,S)~ 4 


L sym Li/2 \ f A R(a) \ 2 2 ] 

h/iKsaJ \a(A,8)J \ j 


(79) 


the limit 5 = 0, AR ^ 0, giving: 


v iv,S =0 

"surf 


\ 


1- 


187 , 


^ AR 2 (ao) 
lP^ Jl/2 aor sat {0) 


A 2 / 3 . (81) 


1/2 


This expression shows that an isovector surface energy can be 
induced in asymmetric nuclei even if no asymmetry is present 
in the bulk. Of course in realistic situations the bulk asymme¬ 
try and the difference between neutron and proton radii are not 
independent variables; in particular the skin is negligeable if 
5 = 0 as we have already assumed in order to obtain Eq. (1781) 
above. 

Eq. ( 1791 ) shows than even in our rather crude approxima¬ 
tion the surface symmetry energy presents a very complex 
dependence on the physical quantities that measure isospin 
inhomogeneity, namely the bulk asymmetry 5 and the neu¬ 
tron skin thickness AR. In particular we find that E™ (A, 5) 
is not quadratic with 5 but has non-negligible linear compo¬ 
nents (see also Fig. [16] below). We have also quantitatively 
tested that both linear and quadratic terms in AR are required 
to correctly reproduce the surface isovector energy. It is in¬ 
teresting to notice that the linear components mix 5 and AR. 
Indeed, as we can see in Eqs. (l80l > and ( 18 1 b . putting to zero 
one of those variables, which both measure the isospin inho¬ 
mogeneities, leads to a quadratic behavior with respect to the 
other variable (cf. eqs (IHOl ) and (l8P). 

Similar to the previous section, the diffuseness is the only 
unconstrained parameter of the model. It can therefore be de¬ 
termined in a variational approach by minimizing the total 
(isoscalar and isovector) surface energy. In section iHBl we 
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FIG. 12. (Color online) Diffuseness as a function of the bulk isospin 
asymmetry. Red lines: Eq. J82b from the minimization of the gaus- 
sian approximation . Blue lines: minimization of the exact numeri¬ 
cally calculated ETF surface energy. Green lines: fit from HF density 
profiles, taken from [20], 


have shown that only the dominant A 2 / 3 terms are important 
to evaluate the diffuseness. For this reason, we neglect again 
terms beyond plane surface, and we approximate the neutron 
skin thickness AR by the hard sphere approximation ARfjs- 
Neglecting the quadratic terms in the expansion in A Rns/a, 
we obtain 


a 2 (A, 8) — aj s (8) 


\ 


n P,„(0) 3J l/2 (S-8*) 


aoARnsiA, 8), 


(82) 


where ajs(S) is the diffuseness obtained in section 1 111 A 11 
by neglecting the isovector component : a/s(5) = 

\/'^rar/(^)/ < ^™r/(^)- We f° un d i n section IIII A II that a/y 
slightly decreases with the isospin asymmetry (see Fig [7]), 
which does not appear consistent with the behavior observed 
in full HF calculations. Now considering in the variational 
principle the isovector term in addition to the isoscalar one, 
the diffuseness a given by Eq. (l82l) acquires an additional term 
which modifies its global 8 dependence. The complete result 
Eq. (|82| > is displayed in Figure [12] We can see that the ad¬ 
ditional term due to the isovector energy contribution inverses 
the trend found section llll A II as expected. More specifically, 
though it does not clearly appear in Eq. ( l82l ). the analytical 
diffuseness is seen to quadratically increase with 5, corrobo¬ 
rating the assumption found in Ref. j2Cl]. 

Although we only considered terms A 2 / 3 , as in a slab 
geometry, the results slightly depend on the nucleus mass as 
shown by the slight dispersion of the different red curves in 
Figure [12] This is due to the neutron skin since AR^gfA. 8 ) 
increases with decreasing mass number A. For comparison, 
the diffusenesses a and a p ^a obtained by a fit of HF density 
profiles in Ref. f20 j ] are also represented in Figure [T2l (green 
curves), as well as the numerically calculated pair (a mm ,a™ m ) 


which minimises the energy (blue curves). 



FIG. 13. (Color online) Fower panel: surface energy per nucleon 
as a function of the bulk isospin asymmetry for four isobaric nucleus 
chains. Full red lines: gaussian approximation using the diffuseness 
Eq. ( 182b . Dash-dotted blue lines: exact numerically calculated ETF 
surface energy using the optimal diffusenesses (a' nm ,a' p m ) (see text). 
Dashed green lines: exact numerically calculated ETF surface energy 
using the diffusenesses from da. 

As we can see, these diffusenesses significantly differ from 
each other, but their consequence on the energy is small as we 
can observe in Fig. [13] which displays the corresponding sur¬ 
face energy E s = E 1 / + E™ per nucleon, for different isobaric 
chains. In this figure, the blue curves correspond to a numer¬ 
ical integration of the ETF energy density, using the diffuse¬ 
nesses which minimize the total surface energy. These results 
can thus be considered as ’’exact” ETF results. The use of 
the very different a and a p values fitted from HF (green lines) 
leads to only slightly different energies, except for the lightest 
isobar chain. The analytical approximation given by the sum 
of Eq. d54l) and Eq. (TT9t . is also plotted (red curves), where 
the diffuseness is given by the analytical formula Eq. (l82l) . 
We can see that our analytical approximation closely follows 
the ’’exact” ETF results. 

All the curves show a positive surface symmetry energy, 
which contrasts with Fig. [TO] As it has been discussed in 1 19T . 
this change of sign is due to the choice between the bulk asym¬ 
metry 8 or the global asymmetry /, in the definition of the bulk 
energy. This choice obviously affects the residual part of the 
energy E s , since the sum of the two gives the same ETF func¬ 
tional. This residual part is, to first order, given by the surface 
symmetry energy as discussed in Ref. |19j]. 

In order to further validate the analytical results of this sec¬ 
tion, quantitative comparisons with Hartree-Fock calculations 
are shown in the next section llll B 31 


3. Comparison to Hartree-Fock calculations 

In this section, we explore the level of accuracy of both 
the no-skin approximation and the gaussian approximation, 
respectively developed in sections UlI B II and IIII B 21 
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As previously discussed, the two different approximations 
lead to two different bulk energetics. Neglecting isospin in¬ 
homogeneities implies that the bulk asymmetry 8 is equalized 
to the average asymmetry I. Thus the bulk quantities p sat and 
E b defined by Eqs. ( l42l ) and (l47l) depend on I, and the total 
energy of a nucleus (A,/) within the no-skin approximation is 
given by 

E N oSkin(A,I) = E b (A,I)+E I s s (A,I)+E™ (A,/), (83) 

where E I S S (A 1 I) is given by Eq. El (with I instead of 5), 
E™ (A,I) by Eq. (l59l >. and the diffuseness is given by Eq. (l63l) . 

On the other hand, the gaussian approximation allows 
defining two independent density profiles. Therefore, the bulk 
energy depends on the bulk asymmetry 8(A,I) defined by 
Eq. ( l46l > and the total energy of a nucleus (A,/) within this 
approximation is given by 

EcaussM = E b (A,8)+E I s s (A,8)+E? (A,8), (84) 

where E ! S S (A,8) is given by Eq. El . E^ V (A,8) by Eq. ( l79l ). 
and the isoscalar diffuseness is given by Eq. ( l8l . 



FIG. 14. (Color online) Total energy E = E b + E s per nucleon as 
a function of the nucleus asymmetry I = 1—2Z/A calculated within 
the no-skin approximation, Eq. 1831 (blue dotted lines) and within 
the gaussian approximation, Eq. 1841 (red full lines), compared to 
nuclear Hartree-Fock energy (stars), a) A = 50; b) A = 100; c) A = 
200; d) A = 400. 

In Figure [14] we compare the analytical expressions ( l8l 
and < [84b with Hartree-Fock energy calculations, for different 
isobaric chains. To compare the same quantities, we used the 
same interaction (SLy4), and we have removed the Coulomb 
energy from the total HF energetics. 

We can see from the figure that the no-skin and the gaussian 
approximations predict close values for the total energy. For 
low asymmetries I < 0.2 where the two models are almost 
undistinguishable: they reproduce the microscopic calcula¬ 
tions with a very good accuracy, especially for medium-heavy 
nuclei A > 100. However, for higher asymmetries I > 0.2 
where the symmetry energy becomes important, a systematic 
difference between the two models appears and increases up 
to ~ 400 keV/A for the highest asymmetries I ~ 0.4: the gaus¬ 
sian approximation is systematically closer to the microscopic 


results than the no-skin model. This observation highlights the 
importance of taking into account the isospin asymmetry in¬ 
homogeneities, considering the neutron skin and at the same 
time differentiating the bulk asymmetry 8 from the global 
one I, as it has been discussed in Ref. [19]. Quantitatively, 
for medium-heavy nuclei, the accuracy of Eq. ( l8l is better 
than ~ 200 keV/A, which is similar to the predictive power of 
spherical Hartree-Fock calculations for this effective interac¬ 
tion, with respect to experimental data. 

To conclude, the gaussian approximation developed in sec¬ 
tion 1 III B 21 provides a reliable analytical formula, especially 
for the surface symmetry energy. For this reason we will only 
use the gaussian approximation to further study the different 
components of the nuclei energetics, as we turn to do in the 
next section. 


IV. STUDY OF THE DIFFERENT ENERGY TERMS 

In this section, we use the analytical formulae based on the 
gaussian approximation detailed in section lHIB 21 to study the 
different components of nuclear energetics. As we have previ¬ 
ously discussed throughout this paper, we can decompose the 
nucleus total energy E into bulk E b and surface E s parts. Both 
can be written as sums of isoscalar Ej s , that is the part inde¬ 
pendent of Pi(r), and isovector Ej v terms. The surface energy 
can be further split into plane surface E sur t °c A 2 / 3 , curvature 
E C urv ^ A 1 3 and mass independent E, m j terms. Finally, we 
can distinguish the local and the non-local /o ,V V/ ' com¬ 
ponents of the surface isoscalar part only, since we did not 
discriminate them in the gaussian approximation used for the 
isovector energy. In summary, the energy of a (A,/) nucleus 
can be written as 


E(A,I)=E b (A,8)+E s (A,8), 

(85) 

E s (A,8)=E I s s (A,8)+E?(A,S), 

(86) 

E 1 / (A, 5) = EZ f (A, 8)+EZ(A, 8 ), 

(87) 

(A, S) = £f r/ (A, 8) + E I c s urv (A, 8)+E I i s nd 

(A, 5), (88) 

Ell f {A,8) = E^ f (A, 8)+E I s s ur N f L (A, 8 ), 

(89) 

Ectrv(A,8) = E^(A,8) +E i c s u % l (A, 8), 

(90) 


where the bijective relation (for a given mass) between I and 
8 is given by Eq. (l46t . The different isoscalar terms £. are 
defined by Eqs. (El to (1571) . with the diffuseness £z(A,<5) de¬ 
termined within the gaussian approximation, Eq. d82i l. The 
isovector components E] v are introduced in Eq. (f79b , where 
the curvature term, in this gaussian approximation, is identi¬ 
cally zero by construction. 

In the following, we will study each of these terms, and 
specifically their dependence with the asymmetry 5. For this 
comparison, we have chosen a representative isobaric chain 
A = 100 for which the ETF approximation was successfully 
compared to HF results in Fig. 1 1 41 for the SLy4 interaction. 
For this choice of mass, 8 ss 0 corresponds to the proton 
dripline and 8 « 0.3 the neutron dripline (see Fig.H. 
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Due to our limited experimental knowledge of the isovector 
properties of the effective interaction, the behavior of the dif¬ 
ferent energy terms with asymmetry is to some extent model 
dependent. In order to sort out general trends we have consid¬ 
ered different Skyrme functionals which approximately span 
the current uncertainties on the density dependence of the 
symmetry energy. 

The corresponding bulk parameters are reported in Ta¬ 
ble U] In this table, the calculated surface coefficients 
(J\/ 2 ,L\/ 2 ,K{/ 2 ) entering Eq. (179b and (l82l) are also given. As 
it is well known |55;], the different interactions are very close 
at half saturation density, reflecting the fact that all Skyrme pa¬ 
rameters have been fitted on ground state properties of finite 
nuclei, which correspond to an average density of the order 
of Piaf/2. Nevertheless, a considerable spread is already seen 
at saturation density, showing that the extrapolation of isovec¬ 
tor properties to unexplored density domains is still not well 
controlled li55ll . 

Concerning the LNS interaction, the parametrization pro¬ 
posed in Ref. corresponds to a too high saturation den¬ 
sity which is not realistic. This induces a trivial deviation 
with respect to the other interactions in both the bulk and sur¬ 
face isovector components. For this reason, only the isovector 
properties of this functional are of interest for this study. 

A more complete study of the effective interactions param¬ 
eter space would be necessary to reach sound conclusions on 
the quantitative model dependence, but from the representa¬ 
tive chosen interactions, we can already dress some qualitative 
interpretations. 

The bulk energy per nucleon is shown in the upper panel of 
Fig-El At low asymmetries, the curves are indistinguishable 
reflecting the good present knowledge of symmetric nuclear 
matter properties. The only exception is given by LNS, which 
presents a global shift with respect to the other functionals. As 
already remarked, this is due to the irrealistically high satura¬ 
tion density of this parametrization (tab. [T]i . However, we can 
see that the behavior with isospin is comparable to the one of 
the other functionals, reflecting a compatible bulk symmetry 
energy. For the highest asymmetries 5 > 0.25, we can see 
that all the parametrizations differ, which reflects the larger 
uncertainties for asymmetric matter. 

The lower panel of Figure Q2] displays the surface correc¬ 
tions. We can see that the qualitative behaviour of the differ¬ 
ent models is the same: E s /A increases with the asymmetry, 
leading to a positive sign of the corresponding symmetry en¬ 
ergy. As it has been already discussed in Ref. 11911 . this comes 
from the consideration of the bulk asymmetry 8 instead of the 
global one I in the definition of the nuclear bulk. 

The increase rate with isospin is not the same in the differ¬ 
ent models, reflecting the different surface symmetry energies 
of the functionals. In particular, the steep behaviour predicted 
by the SkI3 parametrization is due to the stiff isovector prop¬ 
erties of this effective interaction (see L sym and K sym in tab.Q}. 
which lay close to the higher border of the presently accepted 
values for these parameters 155.1. 

Moreover, the four considered interactions predict very dif¬ 
ferent values of E s . In particular, at 5 = 0 for which the SLy4, 
SkI3 and SGI models are in perfect agreement on the bulk en- 
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FIG. 15. (Color online) Bulk (upper panel) and surface (lower 
panel) energy per nucleon as a function of the bulk asymmetry 8 
for isobaric nuclei A = 100, predicted by Eq. ( 184b . Different Skyrme 
interactions are considered: SLy4 l23ll ( full red), SkI3 0 (dashed 
green), SGI IbTIl (dotted blue), LNS l62ll (dashed-dotted black). 


ergy, they however differ from ~ 500 keV per nucleon on the 
surface energies. We will come back to this surprising result 
later in this section. 

Fig- E] shows the energy decomposition of Eqs. < 1 86 b and 
©. As expected, at 8 = 0, though not identically zero (see 
Eq. (ED), the isovector energy (lower panel) is completely 
negligible. This a-posteriori justifies the assumption Ej v (0) = 
0 we made in order to obtain ciq in Eq. ( 1771 ). However, for 
asymmetric systems, though smaller than the isoscalar energy 
(upper panel), the isovector energy cannot be neglected. In¬ 
deed, its dependence with 8 is much stronger, meaning that 
the isovector term is the most important term determining the 
surface symmetry energy . Concerning the mass independent 
term, we can see that it is negligible compared to the other 
components, as expected for the medium-heavy nucleus con¬ 
cerned by this picture. Finally, we can observe that the isovec¬ 
tor energy is not quadratic with 8, thus confirming that the 
linear terms of Eq. < f79b cannot be neglected. 

Fig. [17] shows the predictions of the different functionals 
concerning the parameters associated to the density profiles, 
namely the diffuseness (upper panel), the neutron skin (mid¬ 
dle panel) and their ratio (lower panel). We can see that, for a 
given asymmetry 8, the spread of the diffuseness values given 
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TABLE I. Bulk and surface nuclear properties for the different Skyrme interactions examined in this paper. 


Interaction 

Psat (0) 
(fm- 3 ) 

m* /m 

K S at 

(MeV) 

Jsym 

(MeV) 

^sym 

(MeV) 

K-sym 

(MeV) 

h/2 

(MeV) 

L l/2 

(MeV) 

K l/2 

(MeV) 

SLY4 [231 

0.1595 

0.595 

230.0 

32.00 

46.0 

-119.8 

22.13 

38.6 

-74.0 

SkI3 [60] 

0.1577 

0.577 

258.2 

34.83 

100.5 

73.0 

18.85 

46.7 

-25.2 

SGI [611 

0.1544 

0.608 

261.8 

28.33 

63.9 

-52.0 

16.75 

38.4 

-29.7 

LNS [62] 

0.1746 

0.826 

210.8 

33.43 

61.5 

-127.4 

21.10 

44.6 

-56.8 
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FIG. 16. (Color online) Isoscalar (upper panel) and isovector (lower 
panel) surface energy per nucleon as a function of the bulk asymme¬ 
try 8 for isobaric nuclei A = 100, predicted by Eq. GUD. Different 
Skyrme interactions are considered: SLy4 [23] (full red), SkI3 da 
(dashed green), SGI ll6ltl (dotted blue), LNS |[62] (dashed-dotted 
black). 


FIG. 17. (Color online) Diffuseness a (upper panel), neutron skin 
thickness A R (middle panel) and the ratio A R/a (lower panel) as a 
function of the bulk asymmetry 8 for the isobaric chain A = 100, 
predicted within the gaussian approximation (see text). Different 
Skyrme interactions are considered: SLy4 [23j] (full red), SkI3 na 
(dashed green), SGI dill (dotted blue), LNS |62 1 (dashed-dotted 
black). 


by Eq. (l82l ) is very important, reflecting the poor knowledge 
of this quantity. These large uncertainties can be understood 
considering that the diffuseness does not seem to affect the 
energy in a systematic way. In particular, though SkI3 and 
SGI models surprisingly give the same diffuseness, the corre¬ 
sponding surface properties systematically differ. Moreover, 
this similarity of the diffuseness cannot be straightforwardly 
linked to any specific interaction property or parameter (see 
tab.|T|. This reflects again the fact that the diffuseness is a del¬ 
icate balance of all energy components, and is determined by 
very subtle competing and opposite effects. 

The middle part of the Figure shows the obvious correlation 


between A R = R — R p and 5. It is clear from this behavior that 
quadratic terms in the neutron thickness cannot be neglected 
to correctly estimate the symmetry energy (see Eq. (1791)). It is 
interesting to observe that the SGI and LNS models give very 
close results for this quantity, and the same was true for the 
isovector part of the surface energy in Figure fl6l above. 

This comes from the fact, already observed in the literature 
|55ll , that A R is mainly determined by the slope of the sym¬ 
metry energy L sym |55|] which are close in the SGI and LNS 
models. Our work confirms that the neutron thickness can be 
viewed as a measurement of the L parameter. Indeed, A R can 
be well approximated using the equivalent hard spheres radii 
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Rhs(S)’ RhsAV, see Eq. ([75}. This means that A R can be 
seen as a function of the saturation density p S nt(8). In turn, 
the saturation density is given by Eq. (142} which at first or¬ 
der is quadratic in 8 2 with the coefficient L sym /K sat . Since 
K sar is relatively well constrained, we then understand why 
A R is mainly determined by L sym . In particular, the neutron 
skin thickness is predicted to be the same in the two specific 
interactions SGI and LNS. Since the surface isovector energy 
Eq. ([79} at a given bulk asymmetry mainly depends on the 
neutron skin, this also explains why we obtain the same ener¬ 
gies for the two models in Fig. [16] 

This essential role of A R to determine the symmetry energy 
is confirmed observing from Fig.[T6]and[T7]that Skyrme mod¬ 
els which predict thicker neutron skin, that is higher L sym , give 
systematically larger values of the isovector surface energy. 

The lower part of FigurelTTlshows the ratio A R/a as a func¬ 
tion of 8. Though it is the quantity which mainly governs the 
behavior of Eq. d79} . it does not constrain the surface isovector 
energy E IV . Indeed, same A R / a from the functionals SkI3 and 
SGI lead to different energies (Fig. |T6] lower panel), corrobo¬ 
rating the above discussion: only the L parameter, or equiva¬ 
lently the neutron skin thickness A R, is relevant to determine 
the isovector contribution. 

This stresses the importance of the experimental measure¬ 
ment of neutron skin thickness as a key quantity for the knowl¬ 
edge of the density dependence of the symmetry energy [55]. 



FIG. 18. (Color online) Decomposition of the local (left) and non¬ 
local (right) part of the isoscalar surface energy per nucleon, into 
its surface and curvature component as a function of the bulk asym¬ 
metry 8 for the isobaric chain A = 100, as predicted by Eq. [84}. 
Different Skyrme interactions are considered: SLy4 |23f]_ (full red), 
SkI3 j60h (dashed green), SGI Ifilll (dotted blue), LNS ImI (dashed- 
dotted black). 

To conclude, we study in Figure QjUthe decomposition into 
local and non-local terms as predicted by the different func¬ 
tionals. Only the isoscalar part of the surface energy is consid¬ 
ered because these different terms are mixed up in the gaus- 
sian approximation we have employed for the isovector com¬ 
ponent. 

Again, we can see that the qualitative behavior of the differ¬ 
ent Skyrme models is the same for each specific term. We can 


then safely conclude that the non-local curvature component 
Eain L can be neglected for medium-heavy nuclei A > 100, but 
the local curvature energy has to be taken into account since 
it represents for these nuclei 10% to 25% of the total surface 
local energy, depending on the interaction choice and on the 
asymmetry 8. 

Concerning the 8 dependence of the isoscalar surface ener¬ 
gies in Fig. HU we can notice that the local and non-local parts 
have opposite behaviors, leading to the rather flat curves ob¬ 
served in Fig.[T6j upper panel. In section HTCl and UTl A 31 we 
have shown that the exact equality E^h = E^J 1 ' (Eq.(f36b) 
is obtained only if both curvature and isovector terms are ne¬ 
glected in the determination of the diffuseness. However, the 
neglect of isovector terms leads to a wrong dependence with 8 
as shown in Fig. [7] Thus, isovector terms cannot be avoided. 

The results of Figure[l8]clearly show that, once these terms 
are consistently added in the variational procedure (Eq. (l82t ). 
the equality E 1 ^ = e[ S u ^ L is completely violated for asym¬ 
metric systems. Therefore the isoscalar energy strongly de¬ 
pends on the neutron skin thickness, even if it is an indirect 
dependence through the diffuseness. This shows that, though 
the energy can be splitted into different terms, these latter can¬ 
not be decorrelated and have to be treated altogether. 

We have already observed in Figure |T5] that the different 
functionals predict very different surface energy at 8 = 0, 
which might be surprising considering that the symmetric nu¬ 
clear properties are supposed to be well constrained by ex¬ 
perimental data. An obvious interpretation would be that the 
discrepancy comes from the surface properties, that is the non¬ 
local gradient terms and the (poorely constrained) diffuseness 
parameter. However, comparing the different values of the 
predicted diffuseness at 8 = 0 from Figure [17] we can see 
that cisci < a lns = <*Ski 3 < QsLy 4 - This inequality sequence is 
not respected for the surface energy E sur f(8 = 0) in Fig. [T8] 
meaning that the difference of surface energies cannot be as¬ 
cribed to the diffuseness. 

The possible dependence on the couplings of gradient and 
spin-orbit terms is also excluded. Indeed, we can see from 
Figure [18] that at 8 = 0, the isovector part is zero by def¬ 
inition and therefore the equality E^'h = E 1 ^ 1 is verified. 
This means that the total surface energy for symmetric bulk is 
E S urf = 2 E su ’ r p which does not depend on the non-local terms 
of the functional, but only depends on the bulk interaction 
coefficients (p sat (0),Co,Cj,C e ff,a) according to Eqs. (fT7} - 

We can conclude that the differences of the total surface en¬ 
ergies observed for 5 = 0, that is nuclei very close to isospin 
symmetry, in Figure Q2] does not come from the non-local 
properties but are intrinsically linked to the bulk interaction 
coefficients (Co,C 3 ,C e //, a), though the SLy4, SkI3 and SGI 
models correpond to compatible isoscalar equations of state 
(that is: compatible values for the saturation density p mr (0), 
bulk energy Eb(8 = 0), compressibility K sal and effective 
mass). This shows that, at variance with the skin thickness 
A R which is strongly correlated to the isovector equation of 
state, the nuclear surface energy very poorely constrains the 
equation of state, even for symmetric or quasi-symmetric nu- 
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clei. 


V. SUMMARY AND CONCLUSIONS 

In this paper we have addressed the problem of the determi¬ 
nation of an analytical mass formula with coefficients directly 
linked to the different parameters of standard Skyrme func¬ 
tionals, in the extended Thomas-Fermi (ETF) approximation 
at second order in ft. The purpose of this effort is twofold. 
On one side, such a formula is useful for astrophysical ap¬ 
plications where extendend calculations are needed covering 
the whole mass table and using a variety of effective interac¬ 
tion to assess the sensitivity of astrophysical observables to 
the nuclear physics inputs |3l]]. On the other side, analytical 
expressions of the different coefficients of the mass formula in 
terms of the Skyrme couplings allow a better understanding of 
the correlation between these couplings and the different as¬ 
pects of nuclear energetics, for the construction of optimized 
fitting procedures of the functionals. 

The modelling of Fermi density profiles proposed in 
Ref. [|20j l allows an (almost) exact analytical evaluation of 
the isoscalar part of the nuclear energy, naturally leading to 
the appearance in the surface energy of a curvature term and 
a constant term independent of the baryonic number. The 
diffuseness of the density profile is variationally calculated 
within the same formalism, and a simple analytical expres¬ 
sion is given. The relative importance of local and non-local 
terms is studied in detail. Non-local energy components arise 
both from gradient and spin-orbit in the Skyrme functional, 
and from the higher ft terms in the Wigner-Kirkwood expan¬ 
sion of the kinetic energy. We show that in the limit of semi¬ 
infinite matter the isoscalar surface energy is A 2 / 3 and solely 
depends on the local terms. This remarkable property al¬ 
ready observed in Ref. [; 12t] is however violated in finite nuclei 
even if spherical symmetry is assumed, and both components 
contribute in a complex way to the determination of the sur¬ 
face energy. However, the huge dispersion observed on the 
value of the surface tension for symmetric nuclei in modern 
Skyrme functionals is essentially due to the local couplings, 
even if these different functionals correspond to comparable 
saturation properties of symmetric nuclear matter. This find¬ 
ing means that nuclear matter properties are not sufficient to 
pin down surface properties of finite nuclei even in the sym¬ 
metric case. 

The extension to isospin asymmetric nuclei is highly non¬ 
trivial. No exact analytical integration of the ETF functional 
is possible in the presence of isospin inhomogeneities, and ap¬ 
proximations have to be done. We have proposed two different 
approximations for the determination of the surface symmetry 
energy. The first approximation consists in completely negect- 
ing the difference between the neutron and proton radius, that 
is the neutron skin A R. The resulting surface energy shows a 
quadratic dependence on the isospin asymmetry /, and con¬ 
sists of local and non-local plane surface, curvature and mass 
dependent terms which are simple generalizations of the ex¬ 
pressions obtained for symmetric nuclei. Surprisingly, this 
crude approximation reproduces very well numerical Hartree- 


Fock (HF) results for all stable nuclei up to asymmetries of the 
order of I ss 0.2, and leads to a relatively limited overestima¬ 
tion of the order of « 400 KeV/nucleon close to the driplines. 

A better approximation is obtained if isospin inhomo¬ 
geneities are accounted for. To this aim, we have introduced 
a different radius for the neutron and proton distributions, as 
well as an explicit difference between the global asymmetry I 
and the asymmetry in the nuclear bulk 5, due to both Coulomb 
and neutron skin effects. In this more general case, to obtain a 
mass formula we make the assumption that the surface energy 
density is peaked at the nuclear surface, and curvature terms 
can be neglected. A reproduction of HF results within « 200 
KeV/nucleon at the driplines is obtained, and simple expres¬ 
sions are given for the surface energy and the surface diffuse¬ 
ness parameter. In particular we show that both linear and 
quadratic terms in 8 and A R are needed to correctly explain 
the surface term. Moreover, within this analytical mass for¬ 
mula, we show that the neutron skin is essentially determined 
by the slope of the symmetry energy at saturation, thus con¬ 
firming earlier numerical results from different groups [55]. 
Conversely, the surface symmetry energy is shown to be due to 
a complex interplay of all different local and non-local terms 
in the energy functional. This implies that constraints on the 
symmetry energy parameters {J S ym,L S y m ,K sym ) from mass mea¬ 
surements might be model dependent and misleading. As a 
further developement of this work, we plan to extend the mass 
formula to the case of neutron-rich nuclei beyond the dripline 
in equilibrium with a neutron (and possibly proton) gas. Such 
a parametrization will allow including modifications of the 
nuclear surface energy due the presence of continuum states 
in nuclear statistical equilibrium models, currently used for 
different astrophysical applications in supernova and neutron 
star physics [31]. A self-consistent inclusion of pairing effects 
in the local density BCS approximation, using consistent cal¬ 
culations for the mean field and gap equation with the same 
energy functional, is also in progress 163f l. 
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Appendix A: The Skyrme effective interaction 


while the D, coefficients, associated to the isovector part of 
the energy, are given by: 


The Skyrme functional for the energy density is ex¬ 

pressed as ll23l l53ll 


JT(r) =Jf(r) + jr 0 (r) + J%(r)+J% ff (r) + 

+ ^/i«(r)+^so(r)+^ig(r), (Al) 


where the kinetic term, the effective mass term, the zero-range 
term, the density-dependent term, the finite-range term, the 
spin-orbit term and the spin-gradient term are respectively 



Jffeff = C e ffPT + D e ffP 2 T 3 , 

J#0 = CqP~ +D 0 P 3 , 

3 = (C 3 p 2 +£>3P3 2 )P“, 

J%fin = Cf in (Vp )“ + Dfi n (Vp 3 ) 2 , 
tffiso = Cso J ' Vp + D so J 3 • Vp 3 , 
J^sg=C sg i 2 +D sg jl (A2) 


and where we have introduced the local isoscalar and isovec¬ 
tor particle densities, kinetic densities and spin-orbit density 
vectors: 


P(r) =P»(r) + p p (r), 
p 3 (r)=p„(r) —p p (r), 

T(r) = T„(r) + T p (r), 

T 3 (r) = T ;j (r)-Tp(r), 

J(r) = J n (r)+J p (r), 

J 3 (r) = J„(r) — J p (r). (A3) 


The coefficients C, in equations (IA2b . associated with the 
isoscalar contribution, are linear combinations of the tradi¬ 
tional Skyrme parameters f,-, x, and VV’o as follows: 


Co = gfoi 
C3 = ^ f3 ’ 

Ceff = Yg [3fl + f 2 (4x 2 + 5)], 

Cfin = gY [9*1 — ^2(4X2 + 5)], 

Cso = 

Cs g = ^[h(l- 2x l)~ t2 (l+ 2x 2 )] , 


Do = --fo[2x 0 +1], 

D * = ~ 4 % f3 [2*3 + 1], 

D e ff — ^[#2(2c2 + 1)-#i(2xi + 1)], 

Dfin = [3fi(2xi + 1) -tz{2x 2 + 1 )], 

Dso = YWo, 

Dsg = ^ [ft -* 2 ] . (A5) 

The semi-classical development in Ti , so-called Extended 
Thomas-Fermi (ETF), provides expressions for the kinetic 
densities and spin-orbit density vectors, that is at the second 
order (J: 

T 9 (r) = %(r) + T lq (r) + 0(h A ), (A6) 

J?( r ) = Jog(r) + J 2? (r) + 0(h A ). (A7) 

The results of nuclear matter calculations give the zeroth order 
and read: 


%(r) = ^( 37 r) 2 / 3 pp(r) 5 / 3 , (A 8 ) 

Jop(r) = 0. (A9) 

The Wigner-Kirkwood expansion gives the second order of 
the kinetic densities development: 

^ (r) = T iq (r) + % 2 q (r) + t 2 “ (r) , (A 10) 

with 



T ni 

*2q 


T 2q 


1 (Vp,) 2 

36 Pq 

1 Vpc/V/g 

6 fq 

1 / 2m \ 2 

2 WJ 


3 a P„ 


± p *k. 

6 P ‘> fq 



\2 Pq 



(All) 


The second order of the Thomas-Fermi approximation for the 
spin-orbit currents Jiqi?) reads 


2m W q 
Jlq ~ ~ ~2 Pq ~T i 
n Jq 


(A12) 


where we have introduced the effective mass coefficients 
f q ( r) = m/m*( r) with m*( r) the effective masses, and the 
spin-orbit potentials W t/ (r) as follows |54]]: 


fq = 1 + -fir ( c effP ±D ef fP3) , (A13) 

W q = Cso?P ± D S0 V p 3 
+ 2CsgJ ± 2DygJ 3 , 


(A4) 


(A14) 
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where ± stand for neutrons (protons). 

In several Skyrme interactions (such as SLy4, Sill, SGII...), 
the spin-gradient term J^g are neglected. Therefore in the fol¬ 
lowing, we take C sg = D sg = 0, which in particular uncouples 
the equations (lA12b and (IA 1 4b . For more general Skyrme in¬ 
teractions, a full treatment of the spin-gradient terms should 
be implemented ll54ll . 

In symmetric matter, we can set the usual following equal¬ 
ities, at every location r 

2p 9 (r)=p(r), 

2r ? (r) = r(r) = T 0 (r) + ^(r) + T 2 n, (r) + Tf(r), 

2J,(r)=J(r), (A15) 


with 


T ° = 5 


3 (37i 
2 


2 \ 2/3 


,5/3 


T4 = 


1 (Vp) 2 

36 p 


3 A P, 


« L WpVf 1 A f 1 /V/V 
^ =Z-V- + 2P-r-^7P -T > 


6 / 6 r f 12 

2 \ 2 


/ 


w = 


_ 1 (2m V 


2\h 2 J 


C\ 


■Vp\ 


f ) 


2m C S0 Vp 

J = - T p -—. 

/i 2 / 


(A16) 

(A17) 


where we have used W 9 = C v „Vp and where we have intro¬ 
duced the effective mass 

f = fq = ] + K P with K = -pj-Cgff. (A18) 

With these formulae, the energy density given by Eqs. (lA2t 
straightforwardly reads 


^[ P ] = h( P ) + 1 ~-f(4 + x?)+c fin (v P ) 2 

+ V^(Vp) 2 , (A19) 

where we have highlighted the local terms 

l l (p) = —f*o+C 0 p 2 + C 3 p a+2 , (A20) 


and where we have gathered the spin-orbit current (C V ,,J • Vp) 
and kinetic density (^/t|°) terms which lead to the definition 
of the spin-orbit potential V so = — \^C 2 0 . 


„(*) 

k 

fir 

0 

1 

2 


1 

0 

n 2 /6 

0 


5/3 

-0.758981245 

1.517431001 

-2.60168706 


2 

-1 

7T 2 /6 

—n 2 /3 


cc -b 2 

-1.10223102 

1.72183325 

-3.59345480 


3 

-3/2 

1/2 + tt 2 /6 



4 

-11/6 

1 + 7T 2 /6 



5 

-25/12 

35/24 + ^ 2 /6 



6 

-137/60 

45/24 + ^ 2 /6 



7 

-49/20 

203/90+^/6 



8 

-363/140 

469/180+^/6 



9 

-761/280 

29,531/10,080+^2/6 



10 

-7,129/2,520 

6,515/2,016 + ^2/6 



tJA 

TABLE II. Values of the coefficients q/ calculated via the equa¬ 
tions of appendix IB II 

The calculations for y 6 N are analytical; numerical otherwise. For 
the specific T)^ which depends on the value of a, that is of the ef¬ 
fective interaction, we show here the result considering the SLy4 in¬ 
teraction (a = 1 /6). The are given up to the 7 th order in the 
spin-orbit Taylor expansion (see text). 


1. General formulae 


The Fermi function F(r) (Eq. ([6])) to any power can be in¬ 
tegrated in any dimension in using the following general for¬ 
mula J22ll : 


I/n.y — 47T 


~4 K 


m + 1 
with m eN, ysl 4 


/•+<*> 

/ d rF r {r)r m (Bl) 

Jo 

, M m/flA+r 

K-» + i)I L 1, ( 5 ) 


V? = 


(-l^fdJ 1 ^- 1 ^-! 

V 2 Jo (l+e-“) r 


u k , (B2) 


and the binomial coefficient (™) =m\/(k\(m — k)\). The val¬ 
ues of the coefficients that have been used for this work are 
given in table dill ). 

Equation (IB2b is an approximation for which the tiny error 
is ~ exp(— R/a). One can observe that 


„( 0 ) _( 0 )_ 1 . (k) (k)_ k <k-i) 

Py+l Py — y > Py+l Py — (*>0)(B3) 


Appendix B: Integrals of Fermi functions 


2. Expressions of a 3D integral as ID integrals 


We give here the formulae useful to analytically integrate 
Fermi functions to some power. 


In this section we express the difference A Iy y = I 2 y h.y 
as a sum of 1-dimensional integrals. 
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The moments of the difference between two one¬ 
dimensional Fermi functions F(x) = (l +e x/ ") to different 
powers /,y can be integrated as 1221 ] 

J + x k AFy J (x)dx = a k+l (^ < f ) -T 1 ^, (B4) 

with AF y ^ = Fy-F?. 

Making the change of variable x = r — R, we can express 
the 3-dimensional integral h.y = /di\F y (r) as a sum of three 
1-dimensional integrals of moments of Fermi functions F(x): 

I 2 y = 4k f (x + R) 2 F y {x)dx + 4k [ (x + R) 2 F y {x)dx 

J-R J-oo 

r-R 

-4 K ( x + R) 2 F r (x)<Lx , (B5) 

J —oo 

where we have used the Chasles formula to get integrals over 
the entire slab-space. Assuming that the bulk is reached in the 
’’negative” region, that is F^(x< -R) = 1, we can express the 
difference of two Fermi functions to different powers 

Aly.y — 4kR 2 j AFy y(x)dx+ 8kR j xAFy y{x)dx 
+ 4K j x 2 AFy r (x)dx. (B 6 ) 

J —oo 

Because of the previous approximation, we have spuriously 
inserted a bulk part in Eq. (IB 6 I) . but with the very tiny error 

~ ^exp (~YR/a) — exp (—yR/a)^. Computing Eq. (IB 6 I) with 
Eq. (lB4b and expanding the radius parameter R as a series of 
( a/Rfjs ) until the third order according to Eq. ©, we finally 


get at the third order in ( a/Rns)■ 



with r sat = RhsA 1 / 3 = (jKp S at) Let us notice that we 
can also obtain equation (1B7I) using the general formula (IB2I) . 


Appendix C: Analytical expression for the surface energy 

We show in this section how equation (1B7I) allows to obtain 
an analytical formula for the symmetric local E{: and non¬ 
local E^ l surface energy which lead to equations JT5t . ( IT 6 l l 
(O and (l22l) . We also detail the gaussian approximations as 
a function of 1 -dimensional integrals. 

1. The isoscalar local energy 

The surface local energy E'y only depends on the density 
profile p(r) = p sat F(r) through h(p) = £ r c y p r (see Eq. ((5) 
for the values of y, c y ), such that 

E s = J dr|/t(p) - y^p j = £t>pLA/ y j ■ (Cl) 

Computing with the equation (IB7I ) for the y- values of h(p) 
(y = 5/3, 8/3, 2 and (a + 2)), we obtain 



n ( 0 ) + Kn n ( 0 ) 
'>5/3 + K P^" 7 1 8/3 


CoPsc 


^3 Psat Ra+2 


——A 2 / 2 ' 

rsat 


K 2 m 


( 1 ) , ( 1 ) 
ts/i + Kp.a'T," 6m 


+ Clp“i l ( Pa+2 ' 


K~ 


2 K~ ( 0 ) 
— ^3' 


*Psat l ?? 8/3 - 


2k 1 




a 


Fsat 


2 

A 1 / 3 


a+1 


+ —CoPsat + C 2 p sat 


( 2 ) 2 K 2 ( 0 ) 

Tail- ~r%+2 



(C2) 


where the values of rjy are given in table [TT] Using (IB3b 

(A:) (A;) 

which gives a relation between and 7 / 5 / 3 , we & et l° ca l 

1 (k) (k) 

energy E nb as a function of rj// and i] a j 2 only (k = 0,1, 2) 

(Eqs. ©. CH and (EU)). 


2. The symmetric non-local energy 

The non-local energy Ef L is the integration of a quadratic 
function in the density gradient such that we can put it on the 
form £ r c y (Vp)‘p7 - 2 (see Eq. (IT4b for the values of y, c y ). 
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Expressing the Fermi gradient function as follows 

Vp(r) = p sat VF(r) ; VF(r) = 1 (E 2 (r) -F(r)), (C3) 
we can write 


as 1-dimensional integrals in order to obtain the isovector sur¬ 
face energy as a function of the nucleus mass and of the effec¬ 
tive interaction parameters. As for the symmetric energy, we 
make the variable change x = r — Rm'- 


f nl 

As 


,r-2 


[ dr £c r (Vp) 2 p' 

J y 

^LcypLj dr[^(V 7+2 — F r+I ) - (V 7+1 -F 7 ) 


r 

E c rP™r 

r 


A/7+2,7+1 ^7+1,7 


(C4) 


Eg 


= 47r.fi/ 


- 47T.fi/ 


j: 

L 


dx(x + /?M)”exp 


dx(x + f?M) 2 exp 



(C8) 


For 7 > 1, using the recursion relation ( IB3b . we have 

^2 iji - (^ffi - 

Vy+2 - W+t - (Pr 0) +1)» 

^+2 - Pm - ('Pni - = Y(^T) K’ + Pr”) > 

(C5) 


Since we are interested in the surface energy, we assume that 
the gaussian Sf (r) is zero at the center of the nucleus, such that 
the second integral in Eq. (IC8b is negligible with an accuracy 
~ exp (— R\ I /{2a 2 )). Then integrating the gaussian moments 
straightforwardly lead to 


£ G = 2 (2tt ) 3/2 C7.fi/(T^ + ct 2 ), (C9) 


which allows simplifying the expression of E^ L once we have 
computed Eq. ( 1C4I ) with Eq. ( 1B7I ): 


e? l 


E c rP.sV 


1 

7(7+1) 


3 ^a 2 / 3 +6 

r sat 


7 ?7° ) + 1 



+ 6 



- r\f- ] 




(C6) 


Looking at the definition of the non-local energy E^ L Eq. (IT4l) . 
one can see that there are terms /“ 1 = (1 + K'p j 1 . In order 
to have an expression in the form of Eq. (IC6b . we need to make 
a Taylor expansion, such that / _1 = £ i=0 (—l)'(K’p)'. Then 
we can straightforwardly compute the non-local energy with 
Eq. ( 1C6I > (with 7 = 1, 2, i + 2 and i + 3) to obtain Eqs. (IT5] >. 
(Ell and (E3. 


3. The isovector energy 


In this section we develop the 3-dimensional gaussian Sf(r) 
integral 


E g =4k 



poo 

= 47 1 d rr 2 £/ exp 

Jo 


( r-R M ) 2 \ 
2a 2 J 


where we have used the expression of the variance as a 
function of the energy density second derivatives (see sec- 
tion llll Bb . To have Eq as a function of the mass, we just need 
to express the gaussian maximum position Rm as a function 
of A. If we assume Rm = R , it reads, using Eq. ©: 


Eg = 2 (27t) 3 / 2 cr.fi/ r 2 at 



2n 2 

— 




(CIO) 


In the general case, if we define A R = Rm — R, we find addi¬ 
tional terms, especially curvature: 


E G = 2{2nf 2 G^r sat 


A 2 P + 2—A^+^- 

r S at 


r 2 

' sat 


27r 2 
IT 
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2n 2 A R 

3 r sat 



2 A' 1 / 3 + 0 



(Cll) 


(C7) 
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